JP2010184040A

JP2010184040A - 超音波診断装置を用いた生体内物体の動態解析方法

Info

Publication number: JP2010184040A
Application number: JP2009030184A
Authority: JP
Inventors: Kenji Amaya; 賢治天谷; Yasutaka Nishii; 康隆西井
Original assignee: Tokyo Institute of Technology NUC
Current assignee: Tokyo Institute of Technology NUC
Priority date: 2009-02-12
Filing date: 2009-02-12
Publication date: 2010-08-26

Abstract

【解決課題】超音波診断装置から得られるエコー像を用いた生体内物体の動態解析方法を提供すること。
【解決手段】前記超音波エコー像のある時刻のエコーデータから前記対象物の測定境界データを抽出するステップと、予め用意された前記対象物の形状モデルの表面形状データ、及び位置・姿勢パラメータを用いて計算境界データを構築するステップと、前記測定境界データと前記計算境界データとのマッチングを行い、マッチング率Ｒを目的関数、前記位置・姿勢パラメータを設計変数とする最大値探索問題を解いて、前記マッチング率Ｒが最大となる位置・姿勢パラメータを同定するステップとを備え、さらに、前記超音波エコー像の所定時間におけるすべてのエコーデータについて前記各ステップの操作を行って各時間における前記対象物の位置・姿勢パラメータを同定し、得られた各エコーデータの位置・姿勢パラメータと表面形状データとを用いて、所定の三次元コンピュータグラフィクスの手法により対象物の動きをコンピュータのディスプレイ上に表示する。
【選択図】図５Ａ

Description

本発明は、超音波診断装置を用いた生体内物体の動態解析方法に関し、特に、超音波診断装置から得られるエコー像のエコーデータ（静止画）と、予め用意した生体内物体の形状モデル画像とのマッチングを行い、エコーデータ中における形状モデルの位置・姿勢パラメータを同定する処理を各エコーデータに対して行い、生体内物体の動態解析を行う方法に関する。

人工膝関節置換術(Total_Knee_Arthroplasty:ＴＫＡ)は変形性膝関節症や慢性リウマチなどの疾患、あるいは重度の骨折などの外傷によって破壊された膝関節を人工材料で置換する手術法である。ＴＫＡは除痛や歩行能力の再獲得を目的とする手術法として近年ますます多用されてきている。
しかし、ＴＫＡの問題点として人工膝関節の耐用年数や破損材料を取り換える際には再手術を行うため患者への負担が挙げられる（図１）。図１は人工膝関節の説明と破壊例を示すものであり、左側は破損したポリエチレンインサート、右側は人工膝関節の各部品の名称を示している。一般にポリエチレンインサートは術後１０〜２０年で約半数が破損してしまうと言われている。また、従来困難であったＴＫＡ後の深屈曲(術後の深屈曲は130°以上と定義されている。)を人工膝関節のデザイン改良により容易に獲得できれば患者の生活の質(Quality_of_Life:ＱＯＬ)は非常に向上する。

これらの問題を解決するには、生体内における人工膝関節の動態解析が必要である。動態解析が実現すれば人の運動に対するポリエチレンインサートの摩耗解析やポリエチレンインサートと大腿骨コンポーネントとの接触点解析を行えるため、人工膝関節の形状設計にフィードバックできる。形状設計の改良により人工膝関節の耐用年数の向上やＴＫＡ後の深屈曲の実現を期待でき、患者の身体的負担の軽減、再手術費用の軽減、ＱＯＬの向上を実現できる。
人工膝関節の動態解析の方法として、骨や人工膝関節内にマーカを挿入し、２方向からのＸ線撮影像から立体写真計測法を用いて当該マーカの三次元位置を算出する方法やＸ線撮影像の単面投影のみから物体の動態解析を行う方法がある。しかし、これらの方法はマーカを生体内へ侵襲させる必要がある、Ｘ線による放射線被曝(表１参照)が起こる、検査費用が高い、解析精度が必要とされている値に達していないなどの問題がある。
以下の表１は放射線被曝による人体の影響を示すものであり、線量の単位はミリシーベルト(mSv)である。

長尾智晴，安居院猛，長橋宏：遺伝的手法を用いた２値図形のパターンマッチング，電子情報通信学会論文誌，1993年３月、Vol.J76-D-No.3 河野高廣，羽石秀昭，鈴木昌彦，守屋秀繁，森慎一郎，遠藤真広:４次元ＣＴデータを用いた膝関節の動態解析，第２４回日本医用画像工学会大会(JAMITAnnual Meeting 2005)，予稿集(CD-ROM)・A36 山崎隆治，渡邉哲，中島義和，菅本一臣，冨田哲也，前田大助，佐藤嘉伸，吉川秀樹，田村進一：Ｘ線透視画像を用いた人工膝関節の三次元動態解析システムの開発；日本放射線技術學會雜誌(社団法人日本放射線技術学会)，Vol.61、No.1(20050120) pp.79-87

本発明は、上述のような事情に鑑み為されたものであり、マーカを生体内へ侵襲させることなく、また、Ｘ線による放射線被曝もなく、検査費用が比較的安く、かつ、解析精度を向上させるため、超音波診断装置から得られる動画を用いた生体内物体の動態解析方法を提供することを目的とする。

本発明は、超音波診断装置から得られた生体内物体（以下、「対象物」という。）の超音波エコー像の時系列データを用いて、プログラムされたコンピュータにより前記対象物の動態解析を行う方法であって、該方法は、前記超音波エコー像のある時刻のエコーデータから前記対象物の測定境界データを抽出するステップと、予め用意された前記対象物の形状モデルの表面形状データ、及び位置・姿勢パラメータを用いて計算境界データを構築するステップと、前記測定境界データと前記計算境界データとのマッチングを行い、マッチング率Ｒを目的関数、前記位置・姿勢パラメータを設計変数とする最大値探索問題を解いて、前記マッチング率Ｒが最大となる位置・姿勢パラメータを同定するステップとを備え、さらに、前記超音波エコー像の所定時間におけるすべてのエコーデータについて前記各ステップの操作を行って各時間における前記対象物の位置・姿勢パラメータを同定することによって達成される。
また、本発明の上記目的は、得られた各エコーデータの位置・姿勢パラメータと前記表面形状データとを用いて、所定の三次元コンピュータグラフィクスの手法により、前記対象物の動きを前記コンピュータのディスプレイ上に表示することにより前記対象物の動態解析を行うことによって達成される。

本発明に係る超音波診断装置を用いた生体内物体の動態解析方法によれば、超音波エコー像を用いるため、人体へのマーカ侵襲や放射線被曝、検査費用の問題は解決する。また、超音波エコー像は対象とする測定範囲は小さいがＸ線撮影像よりも解像度が高いため、Ｘ線撮影像から動態解析を行う場合よりも高い精度を期待できる。

人工膝関節の説明と破壊例を示すものであり、左側は破損したポリエチレンインサート、右側は人工膝関節の各部品の名称を示している。人工膝関節を模した簡易模型を示すものである。本発明に用いる簡易模型の形状モデルを示すものであり、上段は簡易模型の全景図、下段左は大腿骨コンポーネント（部品１）、下段中はポリエチレンインサート（部品２）、下段右は脛骨ステム（部品３）を模したものである。人工膝関節の各部品の全景図（左）と各メッシュデータ（右）を示すものであり、図上段は大腿骨コンポーネント（部品１）、図中段はポリエチレンインサート（部品２）、図下段は脛骨ステム（部品３）である。本発明に係る方法の手順を示すフローチャートの一例を示すものである。フローチャートを絵で表したものである。ある時刻tのフレーム画像US_im(t)を示すものである。フレーム画像から測定対象領域を抜き出すところを示す図である。抜き出した測定対象領域に対してエッジ抽出を行っている図である。探索空間におけるマッチング率の特徴を示す図である。簡易模型の超音波画像を示す図である。エッジ抽出の方法を示したものである。距離変換されて平滑化された画像の例を示すものである。距離変換する前の断面図を示すものである。距離変換をした後の断面図の画像例を示すものである。距離変換関数（指数関数）の例を示す図である。Ｘ軸方向の変位量（ズレ）とマッチング率との関係を表したグラフである。Ｙ軸方向の変位量（ズレ）とマッチング率との関係を表したグラフである。各部品のメッシュデータと各部品に対して位置・姿勢パラメータを定めたときの例を示すものである。形状モデルのメッシュの断面データ（左図）と、復元された計算輪郭データ（右図）を示す図である。膝のモデルを用いた座標表現を示す図である。超音波画像を用いた場合の解の不安定さを説明するための図である。チコノフの正則化法を用いない場合の部品１のマッチング率の変化を表したグラフである。チコノフの正則化法を用いない場合の部品２及び３のマッチング率の変化を表したグラフである。横軸にλの値の常用対数、縦軸にチコノフの正則化項をとった片対数グラフであり、チコノフの正則化係数λを導出するためのグラフである。チコノフの正則化法を用いた場合の部品１のマッチング率の変化を表したグラフである。チコノフの正則化法を用いた場合の部品２及び３のマッチング率の変化を表したグラフである。動態解析の結果を示すグラフである。３ＤＣＧを用いた動態解析結果を示す図である。

本発明は、超音波診断装置から得られる超音波エコー像のエコーデータ（静止画）と、予め用意した生体内物体の形状モデル画像とのマッチングを行い、エコーデータ中における形状モデルの位置・姿勢パラメータを同定する処理を各エコーデータに対して行い、生体内物体の動態解析を行う方法に関するものであるが、以下、生体内物体が人工膝関節である場合を例として、本発明に係る方法を説明する。
なお、実施例においては、説明の便宜上、２次元の場合について説明する。
３次元超音波診断装置では２次元断層データを重ねることにより体積的なデータを得る。
そのデータのことをボクセルデータと呼ぶが、ボクセルデータは空間をさいの目に細かく区切ったもので，ちょうど２次元画像でいうピクセルを３次元に拡張したものに相当する。
２次元の繰り返し処理をしても３次元の処理は実現可能である。また後述のようなi,jの添え字は２次元ピクセルデータを表現しているが、添え字をi,j,kなどのように３つ設ければ全く同様に３次元のボクセルデータに対して処理をおこなうことができる。
また、対象物の境界は、２次元では「対象物の輪郭」、３次元では「対象物の表面」となるが、請求項１では、これらを総称して「対象物の境界」と表すこととする。以下の実施例において、「画像」とあるのは、前述のエコーデータの２次元の場合の表現である。

動態解析には人工膝関節の表面形状のメッシュデータと人工膝関節の各部品の位置・姿勢パラメータを用いて形状モデルを構築し、そこから得られる計算輪郭データと超音波診断装置から得られる測定データに対し画像処理を行って得られる測定輪郭データを用いた。このとき各部品の位置・姿勢パラメータによる解空間において、２つの輪郭データの重なりを評価値としたときの最大値を与える点を探索する最大値探索問題として考える。この処理を超音波診断装置から得られる動画の各静止画（フレーム画像）に対して行うことで人工膝関節の動態解析を行うことができた。ここでは人工膝関節のモデルを用いて本発明に係る方法の有効性を検証したが、今後、他の生体内物体の表面形状データがあり、超音波画像を撮れる生体内物体であれば、部位を問わず、その動態解析を行うことが可能になる。なお、解析対象の生体内物体の表面形状データは、その生体内物体をＣＴやＭＲＩで予め撮影し三次元解析を行うことによって作ることができる。

また、ここでは人工膝関節そのものを用いずに、図２に示した人工膝関節を模した簡易模型を用いて本発明に係る方法の有効性を検証した。図３は本発明に用いる簡易模型の形状モデルを示すものであり、上段は簡易模型の全景図、下段左は大腿骨コンポーネント（部品１）、下段中はポリエチレンインサート（部品２）、下段右は脛骨ステム（部品３）を模したものである。簡易模型にはφ＝10mmのアルミ球１個、φ＝30mmのアルミ球２個、16×10.5×60mmの凹型のアクリル棒１本、20×26.5×60mmの凹型のアクリル棒１本を使用した。部品１の大腿骨コンポーネントを３つのアルミ球で模し，部品２のポリエチレンインサートを16×10.5×60mmのアクリル棒で模し、部品３の脛骨ステムを20×26.5×60mmのアクリル棒で模した。また、図４は人工膝関節の各部品の全景図（左）と各メッシュデータ（右）を示すものであり、図上段は大腿骨コンポーネント（部品１）、図中段はポリエチレンインサート（部品２）、図下段は脛骨ステム（部品３）である。

図５Ａは本発明に係る方法の手順を示すフローチャートであり、これに基づいて説明する。なお、図５Ｂはフローチャートを絵で表したものである。
［手順１］超音波診断装置を用いて動態解析を行おうとする部位の超音波断面動画を撮る（ステップＳ１）。
［手順２］超音波断面動画の各フレーム画像から測定輪郭データを得る（ステップＳ２）。
(a)超音波断面動画の各フレームに対してエッジ抽出を行う。
(b)各フレームのエッジ抽出データに対して距離変換を行う。
(c)各フレームのエッジ抽出・距離変換後のデータを測定輪郭データとして保持する。
［手順３］形状モデルから計算輪郭データを得る（ステップＳ３）。
(a)動態解析の対象となる部位の各部品のメッシュデータを用意する。
(b)各部品のメッシュデータを位置・姿勢パラメータを用いて絶対座標系に配置する。
(c)超音波診断装置の測定断面を模擬した位置での断面データを得る。
(d)断面データに対して手順２(a)のエッジ抽出と同じ要領でエッジ抽出する。
(e)エッジ抽出したデータを計算輪郭データとして保持する。
［手順４］手順２と手順３で得られた測定輪郭データと計算輪郭データを用いてマッチングを行う（ステップＳ４）。
［手順５］目的関数をマッチング率、設計変数を形状モデルの各部品の位置・姿勢パラメータとする最大値探索問題を遺伝的アルゴリズム（ＧＡ）を用いて解く（ステップＳ５）。
［手順６］上記手順３、４、５の操作を各フレームの測定輪郭データに対して行い、各フレームにおける各部品の位置・姿勢パラメータを求める（ステップＳ６）。
［手順７］横軸をフレーム数・縦軸を各部品の位置・姿勢パラメータのグラフを用いて測定部位の動態解析を行う（ステップＳ７）。

ここでは、手順１で得られた超音波断面動画の各フレーム画像から測定輪郭データを得る方法、すなわち手順２について説明する。
まず、超音波診断装置から得られる動画をUS_mvとし、動画のある時刻tのフレーム画像をUS_im(t)とし、次式（２．１）に示す。
ここで、t₀、…、t_nは各フレーム画像の時刻とする。ある時刻tのフレーム画像US_im(t)を図６に示す。

フレーム画像US_im(t)に対して画像処理を行い測定輪郭データを得る過程を以下に示す。
図７はフレーム画像から測定対象領域を抜き出している図、図８は抜き出した測定対象領域に対してエッジ抽出を行っている図である。図８の測定エッジデータをUS_im,ed(t)とする。
しかし、得られた測定エッジデータおよび図５Ｂ又は図１９の右に示す計算輪郭データは２値画像であるため、マッチングを行う際のマッチング率(二つの画像の一致率)は探索空間において非常に鋭いピークを持つ場合が多いと予想できる(図９)。
このように、マッチング率が鋭いピークを持つことは、最大値探索を困難にしている。なぜなら、最大値を与える点に近い点に探索点があるとき、近くに最大値を与える点が存在することを知ることは困難だからである。言い換えれば、ピークとなる点から少しでもずれると、マッチング率が極端に落ちるため、最大値を与える点を見つけることが困難になるということである。本発明における探索空間では適応度の分布が不連続であり、適応度の高い遺伝子に近い遺伝子をもつ個体が、同様に高い適応度を持っているとは限らない（非特許文献１参照）。このため、図８で得られた測定エッジデータに対して距離変換を施して“ぼかす”ことで計算輪郭データと測定輪郭データが若干ずれていても比較的高い適応度が得られるようにする。エッジ抽出と距離変換の方法については次以降の段落で説明する。

この段落ではエッジ抽出の方法について説明する。図１０より、簡易模型の超音波画像において水は超音波画像では感知できないため、超音波プローブ(図の上端)から簡易模型の表面までの空間は明瞭に映し出されている。これより、本例ではフレーム画像のエッジ抽出には閾値を用いることにした。図１０の各点は８bitで構成されており、０が黒、２５５が白を表している。
図１１はエッジ抽出の方法を示したものである。図のように超音波画像の左の列から右の列までの各列の上から下に向かって各点の値を見ていき最初に閾値（１０）を超えた点をエッジとする。すると、図１１の右側のグレーの点がエッジとして抽出される。

次に、手順２(b)の距離変換の方法について簡単な例を用いて説明し（非特許文献１参照）、距離変換の有効性について説明する。
まず、距離変換の方法について簡単な例を用いて説明する。２値原画像中の白画素の階調値を０、黒画像の階調値をＬ(Ｌは正の整数)とし、画像中の各白画素w_i(wx_i，wy_i)(i=1,2,・・・)に対して次の［１］〜［３］の処理を行う。

［１］wx_i−L+1≦bx_j≦wx_i+L−1かつwy_i−L+1≦by_j≦wy_i+L−1を満たす、w_i(wx_i,wy_i)の近傍領域内の黒画素b_j(bx_j,by_j)(j=1,2,・・・,n_i)を検出する。近傍領域内に黒画素が存在しない場合は白画素w_i(wx_i,wy_i)の処理を終了し、次の白画素について同じ処理を繰り返す。
近傍領域内に黒画素が存在する場合は次の処理を行う。
［２］次式（２．２）で定義される、w_iとその近傍内の黒画素b_jとのチェス盤距離(８連結距離)D_i,jを次式によりそれぞれ求める。
［３］白画素w_i(wx_i,wy_i)の階調値を、Lm−in{D_i,j}に変更する。原画像は、階調値０からＬまでのＬ＋１段階の階調をもつ階調画像に変換される。図１２は距離変換されて平滑化された画像の例を示しており、Ｌ＝２及びＬ＝３の場合を示している。
測定エッジデータUS_im,ed(t)(図１３)に対して距離変換をした時に得られる画像例を図１４に示す。ただし、距離変換には図１５のように100×２^-rで重み付けを行い、階調値を１００とした。この距離変換後の測定データをUS_im,ed,Ld(t)とする。dは階調値とする。

次に、距離変換の有効性について説明する。図１６は正解値からの差ｄｘの値を−25mm〜+25mmまで0.1mmのステップで動かした時にできるマッチング率の変化を表したグラフである。図１７は正解値からの差ｄｙの値を−25mm〜+25mmまで0.1mmのステップで動かした時にできるマッチング率の変化を表したグラフである。両方のグラフのＬは階調値を表し、L=1,10,100の３種類でマッチング率がどのように変化するかを検証した。図１６、図１７の両図からＬが大きくなるにつれて、マッチング率の鋭いピークが鈍ってきているのがわかる。これより、距離変換関数を用いると測定輪郭データと計算輪郭データが若干ずれていても比較的高い適応度が得られるようになる。

次に、形状モデルから計算輪郭データを求める手順３について説明する。
まず、動態解析の対象となる部位の各部品のメッシュデータを用意する。本実施例における人工膝関節の簡易モデルのメッシュデータはPatranとMathematicaを用いて作成したが、実際の人工膝関節の場合は、それを設計した際のＣＡＤデータを用いてもよい。
各部品のメッシュデータをMsh₁,Msh₂,Msh₃とする(図３)。このとき、各部品の位置・姿勢パラメータを(x_n,y_n,z_n,α_n,β_n,γ_n)(n=1,2,3)と定めると、人工膝関節の位置・姿勢は一つに定まる。これを、式で表す。まず、ｘ軸を中心とした回転を表す関数をR¹(α)、ｙ軸を中心とした回転を表す関数をR²(β)、ｚ軸を中心とした回転を表す関数をR³(γ)とすると下記のように表せる。
これを用いると人工膝関節の位置・姿勢は下記の式（２．６）に表せる。
ここで、
は一つに定まったメッシュデータの各ノードの絶対座標、
はメッシュデータの中心から各ノードへの相対位置、
はメッシュデータの中心の位置の絶対座標、α_n、β_n、γ_nはそれぞれｘ、ｙ、ｚ軸周りの回転角度である。図１８に各部品のメッシュデータと各部品に対して位置・姿勢パラメータを定めたものを示す。

ここでは、計算輪郭データへの復元方法を説明する。計算輪郭データの断面は配置した絶対座標系のｘｙ、ｚｙ平面に平行な平面を用いる。図１９に計算輪郭データへの復元方法の様子を示す。図１９の左図はz=+15mmの断面と断面を横切るメッシュの交点をプロットしたもの、右図は左図の各交点間を線形補完して、前述と同様にエッジ抽出を行い、計算輪郭データを構築したものである。この計算輪郭データをS_im,ed(t)とする。

次に、手順４においてマッチング率を計算する。
図１４よりUS_im,ed,Ld(t)はm×n行列として扱える。mは図１４の縦のpixel数、nは図１４の横のpixel数を表している。以後、US_im,ed,Ld(t)はUと表記する。また、図１９よりS_im,ed(t)は点列データなのでl（小文字のエル）×2の行列として扱える。l（エル）は横軸のpixel数を表す。以後、S_im,ed(t)はSと表記する。
これを式で表すと以下のようになる。
ここで、i=1,2,…,mは図１４の縦軸の大きさ、j=1,2,…,nは図１４の横軸の大きさ、f(i,j)は図１４の横ｉ、縦ｊの位置の階調の大きさ、k=1,2,…,lは図１９の横軸の大きさ、g(k)は図１９の縦軸の大きさを表している。以上より、マッチング率Ｒは下記のように表せる。

次に、目的関数をマッチング率、設計変数を形状モデルの各部品の位置・姿勢パラメータとする最大値探索問題を遺伝的アルゴリズム（ＧＡ）を用いて解く（ステップＳ５）。制約条件は下記に示すとおりである。
＜目的関数＞
＜設計変数＞
＜制約条件＞
部品同士は緩衝しない
人工膝関節の可動域
超音波診断装置のプローブを当てる位置
超音波診断装置のプローブを当てる数
部品２と３は一つの剛体として扱う etc・・・

上記の条件下でＧＡを用いてマッチング率Ｒが最大値をとる各部品の位置・姿勢パラメータを同定する（ステップＳ６）。ここで、得られるパラメータはある時刻tにおける値なので、動態解析を行うためにはt₀・・・t_nの各々について位置・姿勢パラメータを同定する必要がある。

＜制約条件＞
最大値探索問題を容易にするために制約条件を導入する。始めに示す通りの制約条件が存在する。本実施例では記述してある制約条件をすべて用いることとする。
（１）部品の数を減らす
人工膝関節は３つの部品に分かれている（図４）。しかし、人工膝関節はポリエチレンインサート（部品２）と脛骨ステム（部品３）が固定されているものと固定されていないものの２種類存在する。ここでは簡単のため部品２と部品３が固定されているものとして解析した。すると、設計変数は１８個から１２個へと減った。また、部品２と３のメッシュデータの中心の位置は部品２の中心を用いた。
（２）超音波診断装置のプローブの位置と数
本実施例では解析によって超音波プローブの当てる位置・数を変えた。位置・数の値は実施者の経験によって定めた。
（３）人工膝関節の可動域
本実施例では人工膝関節の可動域に制限を設けた。膝の運動の中で一番大きいのは膝の屈伸運動である。そのため、屈伸運動を一番表現しやすい膝のモデルを考えると図２０のように仮想関節中心Ｊを定め，そこから各部品のx,y座標を極座標系で表現する。すると膝の屈伸運動を表現する時にｒは小さく変化し、θは大きく変化する。つまり、極座標系を導入することで探索する際にθに注力することで解を見つけやすくなる。これより、x₁,y₁,x₂,y₂を式で表すと下記のようになる。
各パラメータの制約条件を表２に示す。
また、簡易模型の初期条件を考慮した各パラメータの制約条件を表３に示す。

＜距離変換関数の決定＞
距離変換の有効性については前述した。ここでは本発明に係る方法で有効な距離変換関数を求める。試した距離変換関数は２種類で一つは式(2.2)の線形関数、一つは図１５の指数関数である。上記の線形関数・指数関数に対してそれぞれの階調を1,10,50,100と設定した。これに基づいて各距離変換関数のパラメータ変化量とマッチング率の変化を調べ、マッチング率のピークが鋭い順から順位を付けて並べたものを表４に示す。
表４より以下のことが判明した。
・線形関数・階調値１と指数関数・階調値１はピークが鋭い。
・線形関数・階調値10と指数関数・階調値10は部品２の平行移動(x₂,y₂,z₂の変化)の時のピークが鋭い。
・線形関数・階調値50と線形関数・階調値100は全体的にピークが鈍っている。
・指数関数・階調値50と指数関数・階調値100は全体的にピークは鋭すぎず・鈍すぎず安定していて、この二つのグラフはほぼ一致している。しかし、階調値50の方が計算時間は短い。
以上より、探索問題を解く上でマッチング率のピークの鋭さがすべてのパラメータに対して安定的でかつ計算時間が短縮できる指数関数・階調値50の距離変換関数を用いるのが一番有効であるとわかった。

＜チコノフの正則化法の導入＞
先行技術である従来のＸ線画像をもちいたマッチング方法では部品の全体の投影画像が得られていたので部品の位置・姿勢パラメータを同定する場合には安定に唯一解が得られていた。しかし超音波画像を用いた場合には、超音波を照射した方向のみの表面のみの情報しか得られないので、測定誤差などにより姿勢・位置パラメータの同定が不安定になる（図２１参照）。このように超音波画像をこの問題に適用したときに生じる特徴的な不安定な問題を解決しなければならない。そこで、この問題を解決するために、チコノフの正則化法を導入した。
ここでは式(3.1)の最適化手法に対して「チコノフの正則化法」を適用する有効性について説明する。
図２２及び図２３は、それぞれ、チコノフの正則化法を用いない場合の部品１、部品２及び３のマッチング率の変化を表したグラフであり、横軸が各パラメータの変化量，縦軸がマッチング率の変移を表している。
図２２（Ａ）〜（Ｃ）は、それぞれ、チコノフの正則化法を用いない場合の部品１のx₁、y₁、z₁の変化量とマッチング率の変化との関係を表し、同図（Ｄ）〜（Ｆ）は、それぞれ、チコノフの正則化法を用いない場合の部品１のα_１、β_１、γ_１の変化量とマッチング率の変化との関係を表している。また、図２３（Ａ）〜（Ｃ）は、それぞれ、チコノフの正則化法を用いない場合の部品２，３のx₂、y₂、z₂の変化量とマッチング率の変化との関係を表し、同図（Ｄ）〜（Ｆ）は、それぞれ、チコノフの正則化法を用いない場合の部品２，３のα_２、β_２、γ_２の変化量とマッチング率の変化との関係を表している。
図２２（Ｄ）よりα₁、図２２（Ｅ）よりβ₁、図２２（Ｆ）よりγ₁、図２３（Ｃ）よりz₂、図２３（Ｄ）よりα₂、図２３（Ｅ）よりβ₂のマッチング率のピークが鈍っているのが確認できる。そのため、この６個のパラメータの解を安定的に求めるためにチコノフの正則化法を適用する。すると、式(3.1)は以下のように書き換えられる。
ここで、λはチコノフの正則化係数を表す。次に、λの最適な値を求める。図２４は横軸にλの値の常用対数、縦軸にチコノフの正則化項
をとった片対数グラフであり、チコノフの正則化係数λを導出するためのグラフである。このグラフより以下のことがわかる。
・λ≦0.0001では大幅に減少している。
・0.0001≦λ≦0.001では不安定である。
・0.001≦λでは安定している。
以上のことから、λは0.0001≦λ≦0.001が最適な値であるとわかる。そこで、今回はλ＝0.0005として扱う。

λ＝0.0005を用いたときの各パラメータの変化量とマッチング率の変化を図２５及び図２６に示した。図２５及び図２６は、それぞれ、チコノフの正則化法を用いた場合の部品１、部品２及び３のマッチング率の変化を表したグラフであり、横軸が各パラメータの変化量，縦軸がマッチング率の変移を表している。
図２５（Ａ）〜（Ｃ）は、それぞれ、チコノフの正則化法を用いた場合の部品１のx₁，y₁，z₁の変化量とマッチング率の変化との関係を表し、同図（Ｄ）〜（Ｆ）は、それぞれ、チコノフの正則化法を用いた場合の部品１のα_１，β_１，γ_１の変化量とマッチング率の変化との関係を表している。また、図２６（Ａ）〜（Ｃ）は、それぞれ、チコノフの正則化法を用いた場合の部品２，３のx₂，y₂，z₂の変化量とマッチング率の変化との関係を表し、同図（Ｄ）〜（Ｆ）は、それぞれ、チコノフの正則化法を用いた場合の部品２のα_２，β_２，γ_２の変化量とマッチング率の変化との関係を表している。
この結果、図２５（Ｄ）よりα₁、図２５（Ｅ）よりβ₁、図２５（Ｆ）よりγ₁、図２６（Ｃ）よりz₂、図２６（Ｄ）よりα₂、図２６（Ｅ）よりβ₂のマッチング率のピークが、チコノフの正則化法を用いない場合よりも用いた方が明らかに鋭くなっているのが確認できた。

以上のステップＳ２〜Ｓ６を動態解析を行う範囲のすべてのフレームについて行い、すべてのフレームにおける各部品の位置・姿勢パラメータを同定することにより、動態解析が可能になる。次に動態解析について述べる（ステップＳ７）。
＜動態解析の方法＞
式(2.11)を目的関数、各部品のパラメータを設計変数とした最大値探索問題を各フレーム画像に対して行い、図２７のようなグラフを得る（非特許文献２参照）。図２７の横軸はフレーム数、縦軸は、左側がx,y,zの平行移動[mm]，右側が回転角α、β、γ[度]を表している。また、図２７の結果より、各フレームの位置・姿勢パラメータがわかると図２８のように３ＤＣＧで人工膝関節を表現できる（非特許文献３参照）。すなわち、各部品の表面形状データ（多くの場合は表面のポリゴンメッシュデータ）および各部品の位置・姿勢パラメータが各時間に対して得られたならば、３次元コンピュータグラフィックスの基本的な手法で部品の動きをコンピュータディスプレイ上に俯瞰的に表示することができるからである（参考文献：3D Computer Graphics，Alan H. Watt，出版社:Addison-Wesley出版）。これを用いることで、患者に埋め込まれた人工膝関節を客観的で詳細な評価を行えるため、人工膝関節の可動域制限や不安定な動きなどの不良運動の原因が解明できる。

以上、計測データが２次元画像の場合の処理について説明したが、対象物の断面を複数の方向から高速にスキャンして３次元の立体画像データを取得し、それぞれの断面データのフレームについて本発明に係る方法を適用し、各断面データについての位置・姿勢パラメータを同定し、複数の断面データについてそれらを組み合わせることにより、３次元の超音波診断装置による動態解析にも適用することが可能である。

Claims

超音波診断装置から得られた生体内物体（以下、「対象物」という。）の超音波エコー像の時系列データを用いて、プログラムされたコンピュータにより前記対象物の動態解析を行う方法であって、該方法は、
前記超音波エコー像のある時刻のエコーデータから前記対象物の測定境界データを抽出するステップと、
予め用意された前記対象物の形状モデルの表面形状データ、及び位置・姿勢パラメータを用いて計算境界データを構築するステップと、
前記測定境界データと前記計算境界データとのマッチングを行い、マッチング率Ｒを目的関数、前記位置・姿勢パラメータを設計変数とする最大値探索問題を解いて、前記マッチング率Ｒが最大となる位置・姿勢パラメータを同定するステップとを備え、さらに、
前記超音波エコー像の所定時間におけるすべてのエコーデータについて前記各ステップの操作を行って各時間における前記対象物の位置・姿勢パラメータを同定することにより前記対象物の動態解析を行うことを特徴とする超音波診断装置を用いた生体内物体の動態解析方法。
さらに、前記得られた各エコーデータの位置・姿勢パラメータと前記表面形状データとを用いて、所定の三次元コンピュータグラフィクスの手法により、前記対象物の動きを前記コンピュータのディスプレイ上に表示するステップを備えることにより、前記対象物の動態解析を行うことを特徴とする請求項１に記載の超音波診断装置を用いた生体内物体の動態解析方法。
前記測定境界データを抽出するステップが、
前記フレーム画像に対してエッジ抽出を行うステップと、
前記エッジ抽出されたデータに対して距離変換を行うステップと、
前記距離変換後のエッジ抽出データを前記測定境界データとして保持するステップと、
を含むことを特徴とする請求項１又は２に記載の超音波診断装置を用いた生体内物体の動態解析方法。
前記計算境界データを構築するステップが、
前記対象物の表面形状データを前記位置・姿勢パラメータを用いて絶対座標系に配置するステップと、
前記超音波診断装置の測定断面を模擬した位置での断面データを得るステップと、
前記断面データに対して、前記エッジ抽出と同じ方法でエッジ抽出を行うステップと、
前記エッジ抽出されたデータを前記計算境界データとして保持するステップと、
を含むことを特徴とする請求項３に記載の超音波診断装置を用いた生体内物体の動態解析方法。
前記最大値探索問題の解法に遺伝的アルゴリズムを用いることにより、解の最適化を図ることを特徴とする請求項１乃至４のいずれかに記載の超音波診断装置を用いた生体内物体の動態解析方法。
前記最大値探索問題の解の安定化を図るために、チコノフの正則化法を用いたことを特徴とする請求項１乃至５のいずれかに記載の超音波診断装置を用いた生体内物体の動態解析方法。
前記対象物が人工膝関節である請求項１乃至６のいずれかに記載の超音波診断装置を用いた生体内物体の動態解析方法。