(57)【要約】【解決課題】超音波診断装置から得られるエコー像を用いた生体内物体の動態解析方法を提供すること。【解決手段】前記超音波エコー像のある時刻のエコーデータから前記対象物の測定境界データを抽出するステップと、予め用意された前記対象物の形状モデルの表面形状データ、及び位置・姿勢パラメータを用いて計算境界データを構築するステップと、前記測定境界データと前記計算境界データとのマッチングを行い、マッチング率Ｒを目的関数、前記位置・姿勢パラメータを設計変数とする最大値探索問題を解いて、前記マッチング率Ｒが最大となる位置・姿勢パラメータを同定するステップとを備え、さらに、前記超音波エコー像の所定時間におけるすべてのエコーデータについて前記各ステップの操作を行って各時間における前記対象物の位置・姿勢パラメータを同定し、得られた各エコーデータの位置・姿勢パラメータと表面形状データとを用いて、所定の三次元コンピュータグラフィクスの手法により対象物の動きをコンピュータのディスプレイ上に表示する。【選択図】図５Ａ
2010184040.tif 000002
【特許請求の範囲】
【請求項１】
超音波診断装置から得られた生体内物体（以下、「対象物」という。）の超音波エコー像の時系列データを用いて、プログラムされたコンピュータにより前記対象物の動態解析を行う方法であって、該方法は、
前記超音波エコー像のある時刻のエコーデータから前記対象物の測定境界データを抽出するステップと、
予め用意された前記対象物の形状モデルの表面形状データ、及び位置・姿勢パラメータを用いて計算境界データを構築するステップと、
前記測定境界データと前記計算境界データとのマッチングを行い、マッチング率Ｒを目的関数、前記位置・姿勢パラメータを設計変数とする最大値探索問題を解いて、前記マッチング率Ｒが最大となる位置・姿勢パラメータを同定するステップとを備え、さらに、
前記超音波エコー像の所定時間におけるすべてのエコーデータについて前記各ステップの操作を行って各時間における前記対象物の位置・姿勢パラメータを同定することにより前記対象物の動態解析を行うことを特徴とする超音波診断装置を用いた生体内物体の動態解析方法。
【請求項２】
さらに、前記得られた各エコーデータの位置・姿勢パラメータと前記表面形状データとを用いて、所定の三次元コンピュータグラフィクスの手法により、前記対象物の動きを前記コンピュータのディスプレイ上に表示するステップを備えることにより、前記対象物の動態解析を行うことを特徴とする請求項１に記載の超音波診断装置を用いた生体内物体の動態解析方法。
【請求項３】
前記測定境界データを抽出するステップが、
前記フレーム画像に対してエッジ抽出を行うステップと、
前記エッジ抽出されたデータに対して距離変換を行うステップと、
前記距離変換後のエッジ抽出データを前記測定境界データとして保持するステップと、
を含むことを特徴とする請求項１又は２に記載の超音波診断装置を用いた生体内物体の動態解析方法。
【請求項４】
前記計算境界データを構築するステップが、
前記対象物の表面形状データを前記位置・姿勢パラメータを用いて絶対座標系に配置するステップと、
前記超音波診断装置の測定断面を模擬した位置での断面データを得るステップと、
前記断面データに対して、前記エッジ抽出と同じ方法でエッジ抽出を行うステップと、
前記エッジ抽出されたデータを前記計算境界データとして保持するステップと、
を含むことを特徴とする請求項３に記載の超音波診断装置を用いた生体内物体の動態解析方法。
【請求項５】
前記最大値探索問題の解法に遺伝的アルゴリズムを用いることにより、解の最適化を図ることを特徴とする請求項１乃至４のいずれかに記載の超音波診断装置を用いた生体内物体の動態解析方法。
【請求項６】
前記最大値探索問題の解の安定化を図るために、チコノフの正則化法を用いたことを特徴とする請求項１乃至５のいずれかに記載の超音波診断装置を用いた生体内物体の動態解析方法。
【請求項７】
前記対象物が人工膝関節である請求項１乃至６のいずれかに記載の超音波診断装置を用いた生体内物体の動態解析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、超音波診断装置を用いた生体内物体の動態解析方法に関し、特に、超音波診断装置から得られるエコー像のエコーデータ（静止画）と、予め用意した生体内物体の形状モデル画像とのマッチングを行い、エコーデータ中における形状モデルの位置・姿勢パラメータを同定する処理を各エコーデータに対して行い、生体内物体の動態解析を行う方法に関する。
【背景技術】
【０００２】
人工膝関節置換術(Total_Knee_Arthroplasty:ＴＫＡ)は変形性膝関節症や慢性リウマチなどの疾患、あるいは重度の骨折などの外傷によって破壊された膝関節を人工材料で置換する手術法である。ＴＫＡは除痛や歩行能力の再獲得を目的とする手術法として近年ますます多用されてきている。
しかし、ＴＫＡの問題点として人工膝関節の耐用年数や破損材料を取り換える際には再手術を行うため患者への負担が挙げられる（図１）。図１は人工膝関節の説明と破壊例を示すものであり、左側は破損したポリエチレンインサート、右側は人工膝関節の各部品の名称を示している。一般にポリエチレンインサートは術後１０〜２０年で約半数が破損してしまうと言われている。また、従来困難であったＴＫＡ後の深屈曲(術後の深屈曲は130°以上と定義されている。)を人工膝関節のデザイン改良により容易に獲得できれば患者の生活の質(Quality_of_Life:ＱＯＬ)は非常に向上する。
【０００３】
これらの問題を解決するには、生体内における人工膝関節の動態解析が必要である。動態解析が実現すれば人の運動に対するポリエチレンインサートの摩耗解析やポリエチレンインサートと大腿骨コンポーネントとの接触点解析を行えるため、人工膝関節の形状設計にフィードバックできる。形状設計の改良により人工膝関節の耐用年数の向上やＴＫＡ後の深屈曲の実現を期待でき、患者の身体的負担の軽減、再手術費用の軽減、ＱＯＬの向上を実現できる。
人工膝関節の動態解析の方法として、骨や人工膝関節内にマーカを挿入し、２方向からのＸ線撮影像から立体写真計測法を用いて当該マーカの三次元位置を算出する方法やＸ線撮影像の単面投影のみから物体の動態解析を行う方法がある。しかし、これらの方法はマーカを生体内へ侵襲させる必要がある、Ｘ線による放射線被曝(表１参照)が起こる、検査費用が高い、解析精度が必要とされている値に達していないなどの問題がある。
以下の表１は放射線被曝による人体の影響を示すものであり、線量の単位はミリシーベルト(mSv)である。
【表１】
2010184040.tif 000003
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】長尾智晴，安居院猛，長橋宏：遺伝的手法を用いた２値図形のパターンマッチング，電子情報通信学会論文誌，1993年３月、Vol.J76-D-No.3
【非特許文献２】河野高廣，羽石秀昭，鈴木昌彦，守屋秀繁，森慎一郎，遠藤真広:４次元ＣＴデータを用いた膝関節の動態解析，第２４回日本医用画像工学会大会(JAMITAnnual Meeting 2005)，予稿集(CD-ROM)・A36
【非特許文献３】山崎隆治，渡邉哲，中島義和，菅本一臣，冨田哲也，前田大助，佐藤嘉伸，吉川秀樹，田村進一：Ｘ線透視画像を用いた人工膝関節の三次元動態解析システムの開発；日本放射線技術學會雜誌(社団法人日本放射線技術学会)，Vol.61、No.1(20050120)  pp.79-87
【発明の概要】
【発明が解決しようとする課題】
【０００５】
本発明は、上述のような事情に鑑み為されたものであり、マーカを生体内へ侵襲させることなく、また、Ｘ線による放射線被曝もなく、検査費用が比較的安く、かつ、解析精度を向上させるため、超音波診断装置から得られる動画を用いた生体内物体の動態解析方法を提供することを目的とする。
【課題を解決するための手段】
【０００６】
本発明は、超音波診断装置から得られた生体内物体（以下、「対象物」という。）の超音波エコー像の時系列データを用いて、プログラムされたコンピュータにより前記対象物の動態解析を行う方法であって、該方法は、前記超音波エコー像のある時刻のエコーデータから前記対象物の測定境界データを抽出するステップと、予め用意された前記対象物の形状モデルの表面形状データ、及び位置・姿勢パラメータを用いて計算境界データを構築するステップと、前記測定境界データと前記計算境界データとのマッチングを行い、マッチング率Ｒを目的関数、前記位置・姿勢パラメータを設計変数とする最大値探索問題を解いて、前記マッチング率Ｒが最大となる位置・姿勢パラメータを同定するステップとを備え、さらに、前記超音波エコー像の所定時間におけるすべてのエコーデータについて前記各ステップの操作を行って各時間における前記対象物の位置・姿勢パラメータを同定することによって達成される。
また、本発明の上記目的は、得られた各エコーデータの位置・姿勢パラメータと前記表面形状データとを用いて、所定の三次元コンピュータグラフィクスの手法により、前記対象物の動きを前記コンピュータのディスプレイ上に表示することにより前記対象物の動態解析を行うことによって達成される。
【発明の効果】
【０００７】
本発明に係る超音波診断装置を用いた生体内物体の動態解析方法によれば、超音波エコー像を用いるため、人体へのマーカ侵襲や放射線被曝、検査費用の問題は解決する。また、超音波エコー像は対象とする測定範囲は小さいがＸ線撮影像よりも解像度が高いため、Ｘ線撮影像から動態解析を行う場合よりも高い精度を期待できる。
【図面の簡単な説明】
【０００８】
【図１】人工膝関節の説明と破壊例を示すものであり、左側は破損したポリエチレンインサート、右側は人工膝関節の各部品の名称を示している。
【図２】人工膝関節を模した簡易模型を示すものである。
【図３】本発明に用いる簡易模型の形状モデルを示すものであり、上段は簡易模型の全景図、下段左は大腿骨コンポーネント（部品１）、下段中はポリエチレンインサート（部品２）、下段右は脛骨ステム（部品３）を模したものである。
【図４】人工膝関節の各部品の全景図（左）と各メッシュデータ（右）を示すものであり、図上段は大腿骨コンポーネント（部品１）、図中段はポリエチレンインサート（部品２）、図下段は脛骨ステム（部品３）である。
【図５Ａ】本発明に係る方法の手順を示すフローチャートの一例を示すものである。
【図５Ｂ】フローチャートを絵で表したものである。
【図６】ある時刻tのフレーム画像USim(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)エッジ抽出したデータを計算輪郭データとして保持する。
［手順４］  手順２と手順３で得られた測定輪郭データと計算輪郭データを用いてマッチングを行う（ステップＳ４）。
［手順５］  目的関数をマッチング率、設計変数を形状モデルの各部品の位置・姿勢パラメータとする最大値探索問題を遺伝的アルゴリズム（ＧＡ）を用いて解く（ステップＳ５）。
［手順６］  上記手順３、４、５の操作を各フレームの測定輪郭データに対して行い、各フレームにおける各部品の位置・姿勢パラメータを求める（ステップＳ６）。
［手順７］  横軸をフレーム数・縦軸を各部品の位置・姿勢パラメータのグラフを用いて測定部位の動態解析を行う（ステップＳ７）。
【００１３】
ここでは、手順１で得られた超音波断面動画の各フレーム画像から測定輪郭データを得る方法、すなわち手順２について説明する。
まず、超音波診断装置から得られる動画をUSmvとし、動画のある時刻tのフレーム画像をUSim(t)とし、次式（２．１）に示す。
【数１】
2010184040.tif 000004
ここで、t0、…、tnは各フレーム画像の時刻とする。ある時刻tのフレーム画像USim(t)を図６に示す。
【００１４】
フレーム画像USim(t)に対して画像処理を行い測定輪郭データを得る過程を以下に示す。
図７はフレーム画像から測定対象領域を抜き出している図、図８は抜き出した測定対象領域に対してエッジ抽出を行っている図である。図８の測定エッジデータをUSim,ed(t)とする。
しかし、得られた測定エッジデータおよび図５Ｂ又は図１９の右に示す計算輪郭データは２値画像であるため、マッチングを行う際のマッチング率(二つの画像の一致率)は探索空間において非常に鋭いピークを持つ場合が多いと予想できる(図９)。
このように、マッチング率が鋭いピークを持つことは、最大値探索を困難にしている。なぜなら、最大値を与える点に近い点に探索点があるとき、近くに最大値を与える点が存在することを知ることは困難だからである。言い換えれば、ピークとなる点から少しでもずれると、マッチング率が極端に落ちるため、最大値を与える点を見つけることが困難になるということである。本発明における探索空間では適応度の分布が不連続であり、適応度の高い遺伝子に近い遺伝子をもつ個体が、同様に高い適応度を持っているとは限らない（非特許文献１参照）。このため、図８で得られた測定エッジデータに対して距離変換を施して“ぼかす”ことで計算輪郭データと測定輪郭データが若干ずれていても比較的高い適応度が得られるようにする。エッジ抽出と距離変換の方法については次以降の段落で説明する。
【００１５】
この段落ではエッジ抽出の方法について説明する。図１０より、簡易模型の超音波画像において水は超音波画像では感知できないため、超音波プローブ(図の上端)から簡易模型の表面までの空間は明瞭に映し出されている。これより、本例ではフレーム画像のエッジ抽出には閾値を用いることにした。図１０の各点は８bitで構成されており、０が黒、２５５が白を表している。
図１１はエッジ抽出の方法を示したものである。図のように超音波画像の左の列から右の列までの各列の上から下に向かって各点の値を見ていき最初に閾値（１０）を超えた点をエッジとする。すると、図１１の右側のグレーの点がエッジとして抽出される。
【００１６】
次に、手順２(b)の距離変換の方法について簡単な例を用いて説明し（非特許文献１参照）、距離変換の有効性について説明する。
まず、距離変換の方法について簡単な例を用いて説明する。２値原画像中の白画素の階調値を０、黒画像の階調値をＬ(Ｌは正の整数)とし、画像中の各白画素wi(wxi，wyi)(i=1,2,・・・)に対して次の［１］〜［３］の処理を行う。
【００１７】
［１］wxi−L+1≦bxj≦wxi+L−1かつwyi−L+1≦byj≦wyi+L−1を満たす、wi(wxi,wyi)の近傍領域内の黒画素bj(bxj,byj)(j=1,2,・・・,ni)を検出する。近傍領域内に黒画素が存在しない場合は白画素wi(wxi,wyi)の処理を終了し、次の白画素について同じ処理を繰り返す。
近傍領域内に黒画素が存在する場合は次の処理を行う。
［２］次式（２．２）で定義される、wiとその近傍内の黒画素bjとのチェス盤距離(８連結距離)Di,jを次式によりそれぞれ求める。
【数２】
2010184040.tif 000005
［３］白画素wi(wxi,wyi)の階調値を、Lm−in{Di,j}に変更する。原画像は、階調値０からＬまでのＬ＋１段階の階調をもつ階調画像に変換される。図１２は距離変換されて平滑化された画像の例を示しており、Ｌ＝２及びＬ＝３の場合を示している。
測定エッジデータUSim,ed(t)(図１３)に対して距離変換をした時に得られる画像例を図１４に示す。ただし、距離変換には図１５のように100×２-rで重み付けを行い、階調値を１００とした。この距離変換後の測定データをUSim,ed,Ld(t)とする。dは階調値とする。
【００１８】
次に、距離変換の有効性について説明する。図１６は正解値からの差ｄｘの値を−25mm〜+25mmまで0.1mmのステップで動かした時にできるマッチング率の変化を表したグラフである。図１７は正解値からの差ｄｙの値を−25mm〜+25mmまで0.1mmのステップで動かした時にできるマッチング率の変化を表したグラフである。両方のグラフのＬは階調値を表し、L=1,10,100の３種類でマッチング率がどのように変化するかを検証した。図１６、図１７の両図からＬが大きくなるにつれて、マッチング率の鋭いピークが鈍ってきているのがわかる。これより、距離変換関数を用いると測定輪郭データと計算輪郭データが若干ずれていても比較的高い適応度が得られるようになる。
【００１９】
次に、形状モデルから計算輪郭データを求める手順３について説明する。
まず、動態解析の対象となる部位の各部品のメッシュデータを用意する。本実施例における人工膝関節の簡易モデルのメッシュデータはPatranとMathematicaを用いて作成したが、実際の人工膝関節の場合は、それを設計した際のＣＡＤデータを用いてもよい。
各部品のメッシュデータをMsh1,Msh2,Msh3とする(図３)。このとき、各部品の位置・姿勢パラメータを(xn,yn,zn,αn,βn,γn)(n=1,2,3)と定めると、人工膝関節の位置・姿勢は一つに定まる。これを、式で表す。まず、ｘ軸を中心とした回転を表す関数をR1(α)、ｙ軸を中心とした回転を表す関数をR2(β)、ｚ軸を中心とした回転を表す関数をR3(γ)とすると下記のように表せる。
【数３】
2010184040.tif 000006
これを用いると人工膝関節の位置・姿勢は下記の式（２．６）に表せる。
【数４】
2010184040.tif 000007
ここで、
2010184040.tif 000008
は一つに定まったメッシュデータの各ノードの絶対座標、
2010184040.tif 000009
はメッシュデータの中心から各ノードへの相対位置、
2010184040.tif 000010
はメッシュデータの中心の位置の絶対座標、αn、βn、γnはそれぞれｘ、ｙ、ｚ軸周りの回転角度である。図１８に各部品のメッシュデータと各部品に対して位置・姿勢パラメータを定めたものを示す。
【００２０】
ここでは、計算輪郭データへの復元方法を説明する。計算輪郭データの断面は配置した絶対座標系のｘｙ、ｚｙ平面に平行な平面を用いる。図１９に計算輪郭データへの復元方法の様子を示す。図１９の左図はz=+15mmの断面と断面を横切るメッシュの交点をプロットしたもの、右図は左図の各交点間を線形補完して、前述と同様にエッジ抽出を行い、計算輪郭データを構築したものである。この計算輪郭データをSim,ed(t)とする。
【００２１】
次に、手順４においてマッチング率を計算する。
図１４よりUSim,ed,Ld(t)はm×n行列として扱える。mは図１４の縦のpixel数、nは図１４の横のpixel数を表している。以後、USim,ed,Ld(t)はUと表記する。また、図１９よりSim,ed(t)は点列データなのでl（小文字のエル）×2の行列として扱える。l（エル）は横軸のpixel数を表す。以後、Sim,ed(t)はSと表記する。
これを式で表すと以下のようになる。
【数５】
2010184040.tif 000011
ここで、i=1,2,…,mは図１４の縦軸の大きさ、j=1,2,…,nは図１４の横軸の大きさ、f(i,j)は図１４の横ｉ、縦ｊの位置の階調の大きさ、k=1,2,…,lは図１９の横軸の大きさ、g(k)は図１９の縦軸の大きさを表している。以上より、マッチング率Ｒは下記のように表せる。
【数６】
2010184040.tif 000012
【００２２】
次に、目的関数をマッチング率、設計変数を形状モデルの各部品の位置・姿勢パラメータとする最大値探索問題を遺伝的アルゴリズム（ＧＡ）を用いて解く（ステップＳ５）。制約条件は下記に示すとおりである。
＜目的関数＞
【数７】
2010184040.tif 000013
＜設計変数＞
【数８】
2010184040.tif 000014
＜制約条件＞
部品同士は緩衝しない
人工膝関節の可動域
超音波診断装置のプローブを当てる位置
超音波診断装置のプローブを当てる数
部品２と３は一つの剛体として扱う  etc・・・
上記の条件下でＧＡを用いてマッチング率Ｒが最大値をとる各部品の位置・姿勢パラメータを同定する（ステップＳ６）。ここで、得られるパラメータはある時刻tにおける値なので、動態解析を行うためにはt0・・・tnの各々について位置・姿勢パラメータを同定する必要がある。
【００２３】
＜制約条件＞
最大値探索問題を容易にするために制約条件を導入する。始めに示す通りの制約条件が存在する。本実施例では記述してある制約条件をすべて用いることとする。
（１）部品の数を減らす
人工膝関節は３つの部品に分かれている（図４）。しかし、人工膝関節はポリエチレンインサート（部品２）と脛骨ステム（部品３）が固定されているものと固定されていないものの２種類存在する。ここでは簡単のため部品２と部品３が固定されているものとして解析した。すると、設計変数は１８個から１２個へと減った。また、部品２と３のメッシュデータの中心の位置は部品２の中心を用いた。
（２）超音波診断装置のプローブの位置と数
本実施例では解析によって超音波プローブの当てる位置・数を変えた。位置・数の値は実施者の経験によって定めた。
（３）人工膝関節の可動域
本実施例では人工膝関節の可動域に制限を設けた。膝の運動の中で一番大きいのは膝の屈伸運動である。そのため、屈伸運動を一番表現しやすい膝のモデルを考えると図２０のように仮想関節中心Ｊを定め，そこから各部品のx,y座標を極座標系で表現する。すると膝の屈伸運動を表現する時にｒは小さく変化し、θは大きく変化する。つまり、極座標系を導入することで探索する際にθに注力することで解を見つけやすくなる。これより、x1,y1,x2,y2を式で表すと下記のようになる。
【数９】
2010184040.tif 000015
各パラメータの制約条件を表２に示す。
【表２】
2010184040.tif 000016
また、簡易模型の初期条件を考慮した各パラメータの制約条件を表３に示す。
【表３】
2010184040.tif 000017
【００２４】
＜距離変換関数の決定＞
距離変換の有効性については前述した。ここでは本発明に係る方法で有効な距離変換関数を求める。試した距離変換関数は２種類で一つは式(2.2)の線形関数、一つは図１５の指数関数である。上記の線形関数・指数関数に対してそれぞれの階調を1,10,50,100と設定した。これに基づいて各距離変換関数のパラメータ変化量とマッチング率の変化を調べ、マッチング率のピークが鋭い順から順位を付けて並べたものを表４に示す。
【表４】
2010184040.tif 000018
表４より以下のことが判明した。
・線形関数・階調値１と指数関数・階調値１はピークが鋭い。
・線形関数・階調値10と指数関数・階調値10は部品２の平行移動(x2,y2,z2の変化)の時のピークが鋭い。
・線形関数・階調値50と線形関数・階調値100は全体的にピークが鈍っている。
・指数関数・階調値50と指数関数・階調値100は全体的にピークは鋭すぎず・鈍すぎず安定していて、この二つのグラフはほぼ一致している。しかし、階調値50の方が計算時間は短い。
以上より、探索問題を解く上でマッチング率のピークの鋭さがすべてのパラメータに対して安定的でかつ計算時間が短縮できる指数関数・階調値50の距離変換関数を用いるのが一番有効であるとわかった。
【００２５】
＜チコノフの正則化法の導入＞
先行技術である従来のＸ線画像をもちいたマッチング方法では部品の全体の投影画像が得られていたので部品の位置・姿勢パラメータを同定する場合には安定に唯一解が得られていた。しかし超音波画像を用いた場合には、超音波を照射した方向のみの表面のみの情報しか得られないので、測定誤差などにより姿勢・位置パラメータの同定が不安定になる（図２１参照）。このように超音波画像をこの問題に適用したときに生じる特徴的な不安定な問題を解決しなければならない。そこで、この問題を解決するために、チコノフの正則化法を導入した。
ここでは式(3.1)の最適化手法に対して「チコノフの正則化法」を適用する有効性について説明する。
図２２及び図２３は、それぞれ、チコノフの正則化法を用いない場合の部品１、部品２及び３のマッチング率の変化を表したグラフであり、横軸が各パラメータの変化量，縦軸がマッチング率の変移を表している。
図２２（Ａ）〜（Ｃ）は、それぞれ、チコノフの正則化法を用いない場合の部品１のx1、y1、z1の変化量とマッチング率の変化との関係を表し、同図（Ｄ）〜（Ｆ）は、それぞれ、チコノフの正則化法を用いない場合の部品１のα１、β１、γ１の変化量とマッチング率の変化との関係を表している。また、図２３（Ａ）〜（Ｃ）は、それぞれ、チコノフの正則化法を用いない場合の部品２，３のx2、y2、z2の変化量とマッチング率の変化との関係を表し、同図（Ｄ）〜（Ｆ）は、それぞれ、チコノフの正則化法を用いない場合の部品２，３のα２、β２、γ２の変化量とマッチング率の変化との関係を表している。
図２２（Ｄ）よりα1、図２２（Ｅ）よりβ1、図２２（Ｆ）よりγ1、図２３（Ｃ）よりz2、図２３（Ｄ）よりα2、図２３（Ｅ）よりβ2のマッチング率のピークが鈍っているのが確認できる。そのため、この６個のパラメータの解を安定的に求めるためにチコノフの正則化法を適用する。すると、式(3.1)は以下のように書き換えられる。
【数１０】
2010184040.tif 000019
ここで、λはチコノフの正則化係数を表す。次に、λの最適な値を求める。図２４は横軸にλの値の常用対数、縦軸にチコノフの正則化項
2010184040.tif 000020
をとった片対数グラフであり、チコノフの正則化係数λを導出するためのグラフである。このグラフより以下のことがわかる。
・λ≦0.0001では大幅に減少している。
・0.0001≦λ≦0.001では不安定である。
・0.001≦λでは安定している。
以上のことから、λは0.0001≦λ≦0.001が最適な値であるとわかる。そこで、今回はλ＝0.0005として扱う。
【００２６】
λ＝0.0005を用いたときの各パラメータの変化量とマッチング率の変化を図２５及び図２６に示した。図２５及び図２６は、それぞれ、チコノフの正則化法を用いた場合の部品１、部品２及び３のマッチング率の変化を表したグラフであり、横軸が各パラメータの変化量，縦軸がマッチング率の変移を表している。
図２５（Ａ）〜（Ｃ）は、それぞれ、チコノフの正則化法を用いた場合の部品１のx1，y1，z1の変化量とマッチング率の変化との関係を表し、同図（Ｄ）〜（Ｆ）は、それぞれ、チコノフの正則化法を用いた場合の部品１のα１，β１，γ１の変化量とマッチング率の変化との関係を表している。また、図２６（Ａ）〜（Ｃ）は、それぞれ、チコノフの正則化法を用いた場合の部品２，３のx2，y2，z2の変化量とマッチング率の変化との関係を表し、同図（Ｄ）〜（Ｆ）は、それぞれ、チコノフの正則化法を用いた場合の部品２のα２，β２，γ２の変化量とマッチング率の変化との関係を表している。
この結果、図２５（Ｄ）よりα1、図２５（Ｅ）よりβ1、図２５（Ｆ）よりγ1、図２６（Ｃ）よりz2、図２６（Ｄ）よりα2、図２６（Ｅ）よりβ2のマッチング率のピークが、チコノフの正則化法を用いない場合よりも用いた方が明らかに鋭くなっているのが確認できた。
【００２７】
以上のステップＳ２〜Ｓ６を動態解析を行う範囲のすべてのフレームについて行い、すべてのフレームにおける各部品の位置・姿勢パラメータを同定することにより、動態解析が可能になる。次に動態解析について述べる（ステップＳ７）。
＜動態解析の方法＞
式(2.11)を目的関数、各部品のパラメータを設計変数とした最大値探索問題を各フレーム画像に対して行い、図２７のようなグラフを得る（非特許文献２参照）。図２７の横軸はフレーム数、縦軸は、左側がx,y,zの平行移動[mm]，右側が回転角α、β、γ[度]を表している。また、図２７の結果より、各フレームの位置・姿勢パラメータがわかると図２８のように３ＤＣＧで人工膝関節を表現できる（非特許文献３参照）。すなわち、各部品の表面形状データ（多くの場合は表面のポリゴンメッシュデータ）および各部品の位置・姿勢パラメータが各時間に対して得られたならば、３次元コンピュータグラフィックスの基本的な手法で部品の動きをコンピュータディスプレイ上に俯瞰的に表示することができるからである（参考文献：3D  Computer  Graphics，Alan  H.  Watt，出版社:Addison-Wesley出版）。これを用いることで、患者に埋め込まれた人工膝関節を客観的で詳細な評価を行えるため、人工膝関節の可動域制限や不安定な動きなどの不良運動の原因が解明できる。
【００２８】
以上、計測データが２次元画像の場合の処理について説明したが、対象物の断面を複数の方向から高速にスキャンして３次元の立体画像データを取得し、それぞれの断面データのフレームについて本発明に係る方法を適用し、各断面データについての位置・姿勢パラメータを同定し、複数の断面データについてそれらを組み合わせることにより、３次元の超音波診断装置による動態解析にも適用することが可能である。
【図５Ａ】
2010184040.tif 000021
【図９】
2010184040.tif 000022
【図１２】
2010184040.tif 000023
【図１３】
2010184040.tif 000024
【図１５】
2010184040.tif 000025
【図２０】
2010184040.tif 000026
【図２１】
2010184040.tif 000027
【図２２】
2010184040.tif 000028
【図２３】
2010184040.tif 000029
【図２５】
2010184040.tif 000030
【図２６】
2010184040.tif 000031
【図２７】
2010184040.tif 000032
【図１】
2010184040.tif 000033
【図２】
2010184040.tif 000034
【図３】
2010184040.tif 000035
【図４】
2010184040.tif 000036
【図５Ｂ】
2010184040.tif 000037
【図６】
2010184040.tif 000038
【図７】
2010184040.tif 000039
【図８】
2010184040.tif 000040
【図１０】
2010184040.tif 000041
【図１１】
2010184040.tif 000042
【図１４】
2010184040.tif 000043
【図１６】
2010184040.tif 000044
【図１７】
2010184040.tif 000045
【図１８】
2010184040.tif 000046
【図１９】
2010184040.tif 000047
【図２４】
2010184040.tif 000048
【図２８】
2010184040.tif 000049