"関節領域のMRI3次元スキャン画像データを記憶した画像記憶部を備えたコンピュータが実行するプログラムであって、上記画像記憶部で記憶した画像データの関節骨領域を含む範囲を指定する範囲指定処理、上記指定した範囲内の関節骨領域の輪郭情報を抽出する輪郭抽出処理、同範囲内の関節骨領域を取得する領域拡張処理、及び、関節骨領域の輪郭情報に対する変形を受付ける可変モデル処理により骨領域を抽出する骨領域抽出ステップと、上記骨領域抽出ステップで抽出した骨領域から軟骨領域のエッジ位置群を検出するエッジ検出処理、及び、検出したエッジ位置群の連結から軟骨領域の境界を推定する最適経路検出処理により軟骨領域を抽出する軟骨領域抽出ステップとをコンピュータに実行させることを特徴とするプログラム。上記骨領域抽出ステップの可変モデル処理において、動的輪郭法により関節骨領域の輪郭情報を得ることを特徴とする請求項1記載のプログラム。上記骨領域抽出ステップの輪郭抽出処理において、異方性拡散フィルタを用いて画像上のノイズを低減する画像平滑化処理を実行することを特徴とする請求項1記載のプログラム。上記骨領域抽出ステップの領域拡張処理において、FastMarching法を用いることを特徴とする請求項1記載のプログラム。上記骨領域抽出ステップの領域拡張処理において、動的輪郭法による輪郭曲線とFastMarching法による領域拡張で得られた領域との論理積で関節骨領域を取得することを特徴とする請求項1記載のプログラム。上記軟骨領域抽出ステップの最適経路検出処理において、Dijkstraのアルゴリズム及びFastMarching法のすくなとも一方に基づいて軟骨の表面位置の候補となるエッジ点群から表面位置を推定することを特徴とする請求項1記載のプログラム。関節領域のMRI3次元スキャン画像データを記憶した画像記憶部と、上記画像記憶部で記憶した画像データの関節骨領域を含む範囲を指定する範囲指定処理、上記指定した範囲内の関節骨領域の輪郭情報を抽出する輪郭抽出処理、同範囲内の関節骨領域を取得する領域拡張処理、及び、関節骨領域の輪郭情報に対する変形を受付ける可変モデル処理により骨領域を抽出する骨領域抽出部と、上記骨領域抽出部で抽出した骨領域から軟骨領域のエッジ位置群を検出するエッジ検出処理、及び、検出したエッジ位置群の連結から軟骨領域の境界を推定する最適経路検出処理により軟骨領域を抽出する軟骨領域抽出部とを備えたことを特徴とする画像処理装置。上記特許文献に記載された技術は、軟骨領域の位置と濃度値の知識のファジイルール、半月板の位置と濃度値の知識のファジイルールを用いて最終的に半月板領域の画像データを抽出するようにしたもので、ファジイ推論を実行するための各種メンバーシップ関数の設定の仕方によって抽出精度が大きく左右される。このようなファジイ推論を用いず、できるだけユーザによる入力量を削減しながら、安定した抽出結果を得ることが求められている。本発明は上記のような実情に鑑みてなされたもので、その目的とするところは、ユーザによる習熟度の影響を極力廃して操作入力の量を軽減させながら、安定して関節の軟骨部分の3次元画像を正確に抽出することが可能なプログラム及び画像処理装置を提供することにある。本発明の一態様は、関節領域のMRI3次元スキャン画像データを記憶した画像記憶部を備えたコンピュータが実行するプログラムであって、上記画像記憶部で記憶した画像データの関節骨領域を含む範囲を指定する範囲指定処理、上記指定した範囲内の関節骨領域の輪郭情報を抽出する輪郭抽出処理、同範囲内の関節骨領域を取得する領域拡張処理、及び、関節骨領域の輪郭情報に対する変形を受付ける可変モデル処理により骨領域を抽出する骨領域抽出ステップと、上記骨領域抽出ステップで抽出した骨領域から軟骨領域のエッジ位置群を検出するエッジ検出処理、及び、検出したエッジ位置群の連結から軟骨領域の境界を推定する最適経路検出処理により軟骨領域を抽出する軟骨領域抽出ステップとをコンピュータに実行させることを特徴とする。本発明によれば、ユーザによる習熟度の影響を極力廃して操作入力の量を軽減させながら、安定して関節の軟骨部分の3次元画像を正確に抽出することが可能となる。以下、本発明の一実施形態について図面を参照して説明する。図1は、膝関節軟骨領域抽出プログラムをインストールしたパーソナルコンピュータ以下PC10のハードウェア構成を示す。各種処理制御を司るCPU11とフロントサイドバスFSBを介してチップセットのノースブリッジ12が接続される。このノースブリッジ12は、さらにメモリバスMBを介してメインメモリ13と、またPCI-Expressバスを介してグラフィックコントローラ14及びグラフィックメモリ15と接続される他、チップセットであるサウスブリッジ16とも接続され、主としてこれらの間での入出力制御を実行する。サウスブリッジ16は、キーボード/マウス18、ビデオエンコーダ19、ハードディスク装置HDD20、ネットワークインタフェースI/F21、及びマルチディスクドライブ22と接続され、主としてこれら周辺回路とノースブリッジ12との間の入出力制御を行なう。上記ハードディスク装置20内に、OSオペレーティングシステムと各種のアプリケーションプログラム、各種のデータファイル等に加えて、膝関節置換手術で使用する膝関節軟骨抽出プログラム等が予めインストールされているものとする。なお、上記ビデオエンコーダ19は、与えられたデジタル値の画像信号からアナログ値の画像信号であるRGBビデオ信号を生成して出力し、ここでは図示しないディスプレイ部に送ることで、画像が表示される。また、上記マルチディスクドライブ22は、例えばCDCompactDisc規格、DVDDigitalVersatileDisc規格、及びBDBlue-rayDisc:ブルーレイディクス規格に則った各種光ディスク媒体の再生と記録が可能であり、後述するMRI装置で取得した断層写真等を記録した光ディスク媒体を再生して読出すことで、MRI装置で取得した、患者の膝関節周辺の断層3次元形状データ以下MRI画像データと称するを入力してハードディスク装置20に記録可能とする。MRI画像データとしては、例えばT2階調、16ビット濃淡画像、512×512画素、スライス間隔0.3mm、画素間隔0.3125mmで、1つのデータセット当たり約300枚程度を使用する。なお、これらPC10を構成する個々の要素は、きわめて一般的な周知の技術であるのでその説明は省略するものとする。次に上記実施形態の動作について説明する。本プログラムの実行の前準備として、予め膝関節置換手術を行なう患部の幹部をMRI装置で撮影し、その撮影によって得られたMRI画像データがハードディスク装置20に記憶されているものとする。第1の動作例図2乃至図4は、大腿骨側の膝関節部周辺画像中、骨領域から軟骨領域を抽出するための処理内容を示す。図2はそのメインルーチンであり、CPU11は当初に処理対象となるMRI画像をハードディスク装置20から読出して図示しないディスプレイ画面上に表示させるステップS101。図5は、このときディスプレイ画面で表示される3次元のMRI画像を例示している。同図では、SagittalサジタルビューCoronalコロナビュー及びAxialビューとして、画像中心が同一座標位置となるように互いに直交する3方位からの画像を示している。こうして抽出の対象となる画像を選択して表示させた状態で、軟骨抽出操作を開始するべく、入力ダイアグラムを併せてディスプレイ画面に表示させるステップS102。図6は、この時上記図5で示した画像の横位置に表示される入力ダイアグラムを例示する図である。同図では、タイトル膝軟骨領域抽出下で2つのタブFemur大腿骨Tibia脛骨のうちのFemur側が選択されて、各種入力数値を表示するためのダイアログボックスが表示されている。こうして対象となるMRI画像と入力ダイアログとが表示された状態で、操作開始とともにまず大腿骨側の対象部位を囲む注目領域をユーザの操作によって指定するステップS103。図7は、MRI画像中でユーザの操作によって注目領域を指定した状態を例示する。同図において、Sagittalサジタルビューにおける注目領域A11、Coronalコロナビューにおける注目領域A12、及びAxialビューにおける注目領域A13は3次元空間内で連動するものとして指定されている。また、Coronalコロナビューにおける注目領域A12、及びAxialビューにおける注目領域A13においては、膝関節部の顆間窩の左右の範囲を示すデフォルト位置が点線により、内外側顆のデフォルトの概形入力位置が一点鎖線によりそれぞれ自動的に入力されて表示されている。次いで基準形状の入力位置指定として、大腿骨内外顆の中央付近を指定するステップS104。さらに、図中に点線で示す顆間窩の左右の範囲指定を行なうステップS105。図8は、注目範囲A12,A13内の一点鎖線で示す内側と外側の概形と点線で示す顆間窩の左右の範囲とをユーザの操作により指定した状態を示す。その後、大腿骨側関節部の内側または外側の基準形状と、同じく外側または内側の基本形状とを入力するステップS106,S107。図9及び図10は、大腿骨側関節部の内側図9と外側図10の基本形状をそれぞれディスプレイ画面上でユーザがキーボード/マウス18での操作によりドットを連続してプロットして入力した状態を例示している。これらの基本形状の入力は、全体の処理を迅速に処理すべく補助的にユーザ操作により行なうもので、後の行程でも変形可能であるために非常に大まかな形状であればよく、厳密に骨形状の輪郭位置を正確に指定する必要はない。以上の入力に基づいて、CPU11は大腿骨骨領域の抽出処理を実行するステップS108。図3は、この大腿骨骨領域の抽出処理の詳細な内容を示すサブルーチンである。なお、抽出処理はすべてSagittal面に沿ったスライス画像を中心として行なう。これは、大腿骨遠位端にある軟骨の厚み方向がSagittal面に対して概ね平行で、画像上で骨との境界が現れ易いことに起因する。その当初には、3方向それぞれの画像に対し、注目領域の画像に対して異方性拡散フィルタに通してノイズ成分を減少させた上で、輝度勾配強度画像とCanny法による輪郭画像をそれぞれ生成するステップS121。異方性拡散フィルタは、エッジ保存型の平滑化フィルタで、原画像上のノイズ成分を低減させるとともに、輝度変化の大きい部分ではその境界を跨った平滑化は抑制される。この作用で、組織の分布とその境界をより明確にできる。計算式としては、で表される。拡散係数関数cは調整可能なパラメータを持つ輝度勾配強度の関数であり、勾配の大きい部分での時間的変化量が小さくなる働きをさせる。以下の2つの式が代表的な拡散係数関数である。ここでは指数関数形式のものを用いるものとし、定数Kでエッジに対する感度を調整する。すなわち、上記ステップS121での処理後、CPU11は用意した全ての大腿骨の画像、具体的にはSagittal面に沿ったスライス画像に対する骨領域の抽出処理が終了したか否かを判断するステップS122。ここでまだ全ての画像に対する処理が終了していないと判断すると、次いでCPU11は新たに未処理のMRI画像データを読出し、輪郭部に沿うように動的輪郭法によって基準形状を変形するステップS123。ここで動的輪郭法は、弾性力のある変形可能な閉曲線の内部エネルギとその経路上の画像エネルギの和を最小化させ、対象の形状を得る方法である。画像エネルギと内部エネルギの和を、次式で表す。エネルギを最小化させるということは、画像による力場とその力を受ける曲線の力学的均衡状態を求めることであり、処理内部では力場により時間的に曲線を変形させている。GradientVectorFlow参照:XuCPrinceJLGradientVectorFlow:ANewExternalForceforSnakesIEEEProcConfOn19976671XuCPrinceJLSnakesShapesandGradientVectorFlowIEEETransactionsonImageProcessingMarch19987359369を用いて力場を設定したが、この力場はエッジ情報に基づいて力の方向がエッジに向かうようになっている。上記ステップS123で基準形状の変形を行なった後、FastMarching法による領域抽出を行なうステップS124。ここでは輪郭中央部から開始し、画像輪郭部で進行速度が遅くなるように設定する。上記FastMarching法は、Eikonal方程式である次式F|∇T|=1を解く数値計算法で、Fをある位置での伝搬速度とした場合、Tはその位置への到達時間として表すことができる。この計算法に初期条件の領域を与え、境界を広げていくことで骨領域を抽出する。伝搬速度は画像の輝度勾配強度との関数とし、勾配強度の低い場所で伝わり易く、高い場所で伝わりにくいように定めた。現状は、以下の形式の速度関数、すなわち、F=expを利用している。適当な正規化を施した輝度勾配強度と定数aを与えることで伝搬速度を調整する。T2強調のMRI画像上では軟骨と骨の輝度差が大きいので、FastMarching法による領域がその境界を越えないようにしている。画像の状態によっては境界が曖昧となる部分も生じるため、FastMarching法のパラメータ、エッジや画像輝度の情報も領域の進行条件に加えて、領域の漏れに対処する。この条件を超えて漏れる部分は、後述する動的輪郭法で得られる領域との合成で概ね削られる。計算時間の削減と境界漏れを抑えるために、到達時間が一定の範囲を超えると処理を終了するものとする。この時間設定は、動的輪郭法の結果からそれを囲む矩形範囲を取り出し、その対角距離の約半分に到達する時間として設定する。上記FastMarching法により処理後、上記動的輪郭法とFastMarching法の共通領域を骨であるものとしてメインメモリ13内部のバッファ領域に追加するステップS125。こうして2つの異なる方法で領域を抽出する理由は、各方法で抽出領域が外れ易い部分を互いに補うためである。FastMarching法は、輝度勾配変化の大きい部分では比較的正確に領域を求めることができるが、輝度勾配変化が小さい部分では領域の漏れを生じやすい。これに対して動的輪郭法は、パラメータの調整で境界の輪郭が纏まるようにできるが、明瞭な境界が近くに複数存在する場合、どの境界を同定するのかを安定して判断できない。これは、輝度勾配が大きい部分ではFastMarching法が機能的に働き、小さい部分では動的輪郭法の抑制効果が有効に働くということであり、2つの方法の結果を合成することで互いの方法で外れた部分を除外するものである。こうしてSagittal面に沿ったスライス画像に関する処理を終えた後、再び上記ステップS122からの処理に戻る。以下同様にして、Sagittal面に沿ったスライス画像に関する処理を順次実行し、全ての画像に対する処理が終了すると、ステップS122でそれを判断し、バッファ記憶している、最終的に得られた骨領域の画像に対し、例えばモルフォジー演算によるオープニング処理による平滑化処理によって骨領域表面の凹凸を軽減させるステップS126。すなわち、上記で示した処理は、上述した如くSagittal面に沿ったスライス画像を主体として実行するため、Sagittal面に垂直な方向の領域表面の連続性が低くなり、細かな凹凸やバリ状の余分な領域を生じる。そのため、上記オープニング処理によりそれら精度の落ちる部分を平滑化して低減させるものである。そして、平滑化した骨領域部分を他の領域と区別するべく画像上でマスク化しステップS127、以上でこの図3による骨領域のサブルーチンを終了し、図2のメインルーチンに戻る。図11及び図12は、上記図9及び図10に対応し、上記注目領域A11~A13内でマスク化した部分をハッチングにより示している。図2のメインルーチンでは、上記したステップS108での骨領域の抽出処理に続き、関節液の輝度を指定するステップS109。図13は、この関節液の輝度指定状態でディスプレイ画面に表示される、マスク生成中の処理課程を示す小ウィンドウを例示している。次いでCPU11は、上記抽出した骨領域に基づいて軟骨領域の抽出処理を実行するステップS110。図4は、この大腿骨軟骨領域の抽出処理の詳細な内容を示すサブルーチンである。なお、この抽出処理も全てSagittal面に沿ったスライス画像を中心として行なう。なお、大腿骨遠位端の軟骨は、前方では膝蓋骨と接触する範囲、後方では顆間窩の範囲をのぞいて内外側顆の上方まで分布している。この分布範囲と軟骨表面の境界を検出して軟骨領域を求める。処理はスライス画像の単位で行ない、骨領域の輪郭周辺の輝度分布、輝度勾配、及び局地的な形状などの条件から、軟骨と大腿骨が接触している境界範囲を判断する。その当初には、図3の処理で得た骨領域の輪郭を取得した上でステップS141、まず全てのSagittal面に沿ったスライス画像に対する処理を終えていないことを確認した後ステップS142、あらためて未処理のSagittal面に沿ったスライス画像を1枚取得して、以下に示す条件を満たす点列から形成される、輪郭中心からの輪郭位置を取得するステップS143。分布範囲の端と見なす条件は、輝度勾配強度の分布をもとに、割り出したしきい値を下回る。側面の輝度勾配の脇が大きく骨内部方向に向かっている。輪郭の進行方向が逆になる。顆間窩範囲の場合Sagittal面における輝度勾配の向きが規定値から外れる。これらの条件で骨表面における軟骨分布の端位置候補点の列を求めた後、上記ステップS142からの処理に戻り、全てのSagittal面に沿ったスライス画像に対する処理を終えていないことを確認する。こうして各スライス画像毎に骨表面における軟骨分布の端位置候補点の列を求め、全ての画像に対する同処理が終了すると、上記ステップS142でそれを判断する。次いで、上記スライス画像毎に求めた軟骨分布の端位置候補点を連結するべく、AxialとCoronalの2方向それぞれに多項式による曲線を近似し、軟骨の端位置を推定する。空間上の経路は、次式P≡,Ycoronal,z但し、x:前後、y:上下、z:左右。として表される。この多項式曲線の最近傍にある輪郭位置を軟骨端点とする。骨前方では全てのスライス画像に対して多項式近似を行ない、骨後方では内外顆と顆間窩ステップS105でユーザが入力した顆間窩範囲の3箇所に分け、それぞれ近似させたものを合成する。近似は最小二乗法で行なう。この工程で、1つのスライス画像に対し、前方と後方一対の端点が得られる。Sagittal面に沿ったスライス画像に対する処理を終えていないことを確認した後ステップS144、これら端点間の大腿骨表面から外向きに一定距離内でエッジ点を検出する。このエッジ情報を生成するために、骨輪郭に軟骨の最大厚み分を加味した範囲の画像に、異方性拡散フィルタを適用後、Canny法によるエッジ検出を行なうステップS145。軟骨の最大厚みは、例えば5mmと仮定する。参照:HallFMWyshakGThicknessofArticularCartilageintheNormalKneeTheJournalofBoneandJointSurgery198062408413骨領域を囲む矩形の中心を原点とした極座標で骨表面から放射方向外向きに軟骨の最大厚み以内のエッジ点を収集するステップS146。このエッジ点群を軟骨表面の境界候補とし、最適経路検出を行なって外周軟骨表面の輪郭線を生成する。この検出にはDijkstraのアルゴリズムに基づいた方法を用いる。参照:DijkstraEWAnoteontwoproblemsinconnecxionwithgraphsNumerischeMathematik19591269271この経路の点群は不連続なので、多項式近似で点の得られなかった部分の補間を行なうとともに、経路の平滑化も施す。得た内外周の内側を軟骨領域としてメインメモリ13内に設けたバッファ領域に追加した後ステップS148、次のスライス画像を処理するべく上記ステップS144からの処理に戻る。こうしてスライス画像毎に同様の処理を繰返し実行し、軟骨表面の輪郭線をメインメモリ13のバッファ領域に蓄積する。そして、全てのスライス画像に対する処理が終了すると、ステップS144でそれを判断する。この時点では上記骨領域抽出処理の場合と同様に、Sagittalスライス面に直交する方向にバリ状のノイズが生じ易いので、バッファ記憶している、最終的に得られた軟骨領域の画像に対して、例えばモルフォジー演算によるオープニング処理による平滑化処理によって軟骨領域表面の凹凸を平均化して軽減させるステップS149。そして、平滑化した軟骨領域部分を他の領域と区別するべく、骨領域とは異なる色により画像上でマスク化しステップS150、以上でこの図4による軟骨領域抽出のサブルーチンを終了するとともに、図2のメインルーチンにおいても大腿骨側の軟骨領域の抽出処理を終了する。図14及び図15は、上記図11及び図12に対応し、上記注目領域A11~A13内で軟骨領域をさらにマスク化した部分を骨領域のハッチングとは異なるハッチングにより示している。図16は、大腿骨側の膝関節で抽出された骨領域の3次元モデルと、その骨領域に対して抽出された軟骨領域の同モデルの画像を示す。図16Aは、抽出した骨領域BNの3次元モデルである。図中の点SPは、顆間窩の中央に位置し、人工関節置換手術を行なう場合にリーミング等を行なう基準点である。図16Bは、上記図16Aの内容に加え、抽出した軟骨領域CTの3次元モデルである。このように、ユーザが入力するのは、対象部位を囲む注目領域、内側部と外側部の概形、及び顆間窩の左右の範囲を示す位置の3項目とし、ユーザによる習熟度の影響を極力廃して操作入力の量を軽減させながら、安定して大腿骨側の膝関節の軟骨部分の3次元画像を正確に抽出することが可能となる。第2の動作例図17乃至図20は、脛骨側の膝関節部周辺画像中、骨領域から軟骨領域を抽出するための処理内容を示す。図17はそのメインルーチンであり、CPU11は当初に処理対象となるMRI画像をハードディスク装置20から読出して図示しないディスプレイ画面上に表示させるステップS301。このときディスプレイ画面で表示される3次元のMRI画像は、例えば上記図5で示したものとほぼ同様であるものとする。図5では、SagittalサジタルビューCoronalコロナビュー及びAxialビューとして、画像中心が同一座標位置となるように互いに直交する3方位からの画像を示している。こうして抽出の対象となる画像を選択して表示させた状態で、軟骨抽出操作を開始するべく、入力ダイアグラムを併せてディスプレイ画面に表示させるステップS302。このとき上記図5で示した画像の横位置に表示される入力ダイアグラムについては、概ね上記図6に示したものと同様であるものとする。但し、上記図6ではタイトル膝軟骨領域抽出下で2つのタブFemur大腿骨Tibia脛骨のうちのFemur側が選択されていたが、この脛骨側の処理においてはTibiaのタブが選択された状態で表示される。こうして対象となるMRI画像と入力ダイアログとが表示された状態で、操作開始とともにまず脛骨側の対象部位を囲む注目領域をユーザの操作によって指定するステップS303。図21は、MRI画像中でユーザの操作によって注目領域を指定した状態を例示する。同図において、Sagittalサジタルビューにおける注目領域A31、Coronalコロナビューにおける注目区領域A32、及びAxialビューにおける注目領域A33は3次元空間内で連動するものとして指定されている。次いで基準形状の概形を入力する行程として、まずSagittal面での基準形状とCoronal面での基準形状とを入力するステップS304,S305。図22及び図23は、脛骨側関節部のSagittal面図22とCoronal面図23の基本形状をそれぞれディスプレイ画面上でユーザがキーボード/マウス18での操作によりドットを連続してプロットして入力した状態を例示している。これらの基本形状の入力は、全体の処理を迅速に処理すべく補助的にユーザ操作により行なうもので、後の行程でも変形可能であるために非常に大まかな形状であればよく、厳密に骨形状の輪郭位置を正確に指定する必要はない。以上の入力に基づいて、CPU11は脛骨骨領域の抽出処理を実行するステップS306。図18は、この脛骨骨領域の抽出処理の詳細な内容を示すサブルーチンである。なお、抽出処理はすべてSagittal面及びCoronal面の2つの方向に沿ったスライス画像を中心として行なう。その当初には、3方向それぞれの注目領域の画像に対して異方性拡散フィルタを用いて輪郭部を残しながらノイズ成分を減少させた上でステップS321、Canny法によるエッジ検出を行ない、輪郭画像を生成するステップS322。上記ステップS322での処理後、CPU11は用意した全ての脛骨の画像で、Sagittal面の注目領域内全スライス画像に対する基準形状の変形のための画像処理が終了したか否かを判断するステップS324。ここでSagittal面の全ての画像に対する処理がまだ終了していないと判断すると、次いでCPU11は新たに未処理のMRI画像データを読出し、輪郭部に沿うように動的輪郭法によって基準形状を変形するステップS324。変形された曲線の内側の領域を骨領域の候補として、メインメモリ13内部のバッファ領域に更新記憶する。こうしてステップS323,S324,S325の処理を繰返し実行し、Sagittal面の全ての画像を用いて輪郭部に沿うような動的輪郭法による基準形状の変形に関する画像処理を実行する。そして、Sagittal面の全ての画像を用いた画像処理を終了した時点で、CPU11がステップS323でそれを判断する。その後、CPU11は用意した全ての脛骨の画像で、Coronal面の注目領域内全スライス画像に対する基準形状の変形のための画像処理が終了したか否かを判断するステップS326。ここでCoronal面の全ての画像に対する処理がまだ終了していないと判断すると、次いでCPU11は新たに未処理のMRI画像データを読出し、輪郭部に沿うように動的輪郭法によって基準形状を変形するステップS327。変形された曲線の内側の領域を骨領域の候補として、メインメモリ13内部のバッファ領域に更新記憶する。この時点で、バッファ領域には既にSagittal面の処理で記録された骨領域候補の符号あり、Coronal面の領域候補の符号と比較し、両方が骨領域候補とする画素のみ候補として残される。Sagittal面とCoronal面の2方向を用いるのは、骨の境界が強く出ない水平方向の末端部の領域を補い合うためである。こうしてステップS326,S327,S328の処理を繰返し実行し、Coronal面の全ての画像を用いて輪郭部に沿うような動的輪郭法による基準形状の変形に関する画像処理を実行する。そして、Coronal面の全ての画像を用いた画像処理を終了した時点で、CPU11がステップS326でそれを判断する。その後CPU11は、3次元一次微分フィルタによる輝度勾配強度の算出、及び3次元FastMarching法による領域抽出を行なうステップS329。3次元FastMarching法では、3次元で処理した輝度勾配強度から伝搬速度を算出させるものとし、領域漏れを少なくするために、輝度がしきい値以上の部分とエッジ点によるマスクを用意する。FastMarching法による領域は上記マスク内には入らない。上記FastMarching法により処理後、上記動的輪郭法とFastMarching法の共通領域を骨であるものとしてメインメモリ13内部のバッファ領域の内容を更新記憶するステップS330。こうして2つの異なる方法で領域を抽出することにより、各方法で抽出領域が外れ易い部分を互いに補うことができる。こうして注目領域の3次元の骨領域画像を更新して記憶した後、さらに最終的に得られた骨領域の画像に対し、オープニング処理としての例えばモルフォジー演算による平滑化処理によって骨領域表面の凹凸を軽減させるステップS331。骨表面の細かな凹凸やバリ状の余分な領域を平滑化して低減させた上で、平滑化した骨領域部分を他の領域と区別するべく画像上でマスク化しステップS332、以上でこの図18による骨領域のサブルーチンを終了し、図17のメインルーチンに戻る。図24は、上記図23に対応し、上記注目領域A31~A33内でマスク化した部分をハッチングにより示している。またこのとき、例えばAxial面の注目領域A33における突起PTのように、抽出した結果が不自然であると判断される骨領域に対しては、手動によりマスク部分を変形して修正することを可能としてもよい。図25は、上記図24で表示されたマスク部分に対し、ユーザが明らかに不自然であると思われる突起PT部分を手動操作により削除した状態を示す。こうした手動修正による結果は、骨領域の抽出結果に反映して更新記憶される。その後に図17のメインルーチンでは、CPU11が抽出した骨領域から軟骨の輪郭を自動検出するステップS307。なお、脛骨近位端の軟骨は、大腿骨の内外顆に対応する位置に顆間隆起を挟んで2つに分かれて分布している。図19は、この軟骨領域の輪郭の自動検出の詳細な内容を示すサブルーチンである。その処理当初にCPU11は、抽出した脛骨近位端の骨領域の画像から軟骨の2つの分布範囲を自動検出するステップS341。次いで、抽出した骨領域から骨表面の高さマップを生成するステップS342。生成した高さマップに異方向性拡散フィルタをかけてノイズ成分を除去した後ステップS343、Sobelフィルタを使って勾配ベクトルの強度を算出するステップS344。その後、算出した勾配ベクトルの強度に基づき、軟骨表面の中心から放射方向に勾配強度が0ゼロ以上となる場所を複数点探してスプライン化することで軟骨表面の外枠の初期状態を取得するステップS345。この軟骨表面の外枠に対し、さらに動的輪郭法を用いて特に勾配強度の強い場所の外枠を変形して、高さ変化の大きい構造である軟骨表面に沿わせるステップS346。上記算出した軟骨表面の外枠を用い、次に軟骨領域の外枠の初期状態を取得するステップS347。この初期状態の軟骨領域の外枠に対し、動的輪郭法を用いて変形して軟骨領域に沿わせることで、2つの軟骨領域外枠を得るステップS346。こうして生成した2つの閉曲線の画像を軟骨領域外枠として出力しステップS347、以上で図19のサブルーチンを終了して、図17のメインルーチンに戻る。図26は、このときCPU11がディスプレイ画面で表示する、自動検出した軟骨領域の外枠の画像を例示するものである。上述した如く脛骨近位端の軟骨は、大腿骨の内外顆に対応する位置に顆間隆起を挟んで2つに分かれて分布している。図17のメインルーチンでは、こうして自動検出した軟骨領域の輪郭に対して、ユーザの手動による分布範囲の修正を受付けるステップS308。図27は、この手動による修正を受付けた軟骨領域の輪郭を示す。上記図26で示した内容と比較しても、不自然な凹凸を削除する一方で、必要に応じて領域を広げている部分があることがわかる。その後、あらためて軟骨領域の外枠を指定した画像について左右いずれの足の脛骨に関するものであるのかを指定した後ステップS309、関節液の輝度を指定するステップS310。次いでCPU11は、脛骨軟骨領域の抽出を行なうステップS311。図20は、この脛骨軟骨領域の抽出処理の詳細な内容を示すサブルーチンである。その当初には、上記ステップS308の処理で修正した軟骨領域の輪郭を含む注目領域を生成した上でステップS361、当該注目領域に関する全画像データを取得するステップS362。取得した画像データに対して異方向性拡散フィルタをかけてノイズ成分を除去した後ステップS363、3方向からCanny法でエッジ点群を抽出して合成するステップS364。次いで、骨領域内に属する、骨内部に存在するエッジ点を除去しステップS365、さらにエッジ点の分布を考慮し、骨から一定以上遠方上方に位置するエッジ点を軟骨に属さないものと推定して除去するステップS366。その後、2つに分かれた軟骨領域の内側または外側のSagittal面に関する処理が終了したか否かを判断する。終了していないと判断された場合、軟骨領域の画像について未処理のSagittal面に沿ったスライス画像と軟骨分布域の交点を始終端としてエッジ情報点群を通る最適経路をFastMarching法を用いて検出するステップS369。参照:SethianJAFastMarchingMethodSIAMReview199841199235、DijkstraEWAnoteontwoproblemsinconnecxionwithgraphsNumerischeMathematik19591269271ここで得られた経路を軟骨領域の内側または外側のSagittal面における軟骨表面位置の候補とする。軟骨領域の内側または外側のSagittal面に関する処理が終了したか否かを判断し、終了したと判断された場合、同様の処理を軟骨領域の内側または外側のCoronal面でも行う。得られた経路を軟骨領域の内側または外側のCoronal面における軟骨表面位置の候補とする。軟骨領域の内側または外側のCoronal面に関する処理が終了したか否かを判断し、終了したと判断された場合、続いて、同様の処理を軟骨領域の外側または内側のSagittal面でも行う。得られた経路を軟骨領域の外側または内側のSagittal面における軟骨表面位置の候補とする。軟骨領域の外側または内側のSagittal面に関する処理が終了したか否かを判断し、終了したと判断された場合、同様の処理を軟骨領域の外側または内側のCoronal面でも行う。得られた経路を軟骨領域の外側または内側のCoronal面における軟骨表面位置の候補とする。次にSagittal面とCoronal面の2方向で軟骨表面位置情報を平均化して合成しステップS372、さらに軟骨表面位置情報に残るノイズ成分をSavizky-Golayによる2次元のローパスフィルタで平滑化して低減させるステップS373。参照:SavizkyAGolayMJESmoothingandDifferentiationofdatabysimplifiedleastsquaresproceduresAnalyticalChemistry19643616271639、GorryPAGeneralleastsquaressmoothinganddifferentiationbytheconvolutionmethodAnalyticalChemistry199062570573そして、ノイズ成分を低減した軟骨表面の位置情報に基づき、骨と軟骨表面の間を軟骨領域とし、オープニング処理で表面の凹凸を軽減した後に、軟骨領域部分を他の領域と区別するべく、骨領域とは異なる色により画像上でマスク化しステップS374、以上でこの図20による軟骨領域抽出のサブルーチンを終了するとともに、図17のメインルーチンにおいても脛骨側の軟骨領域の抽出処理を終了する。図28は、上記図25に対応し、さらに軟骨領域をマスク化した部分を骨領域のハッチングとは異なるハッチングにより示している。図29は、脛骨側の膝関節で抽出された骨領域の3次元モデルと、その骨領域に対して抽出された軟骨領域の同モデルの画像を示す。図29Aは、抽出した骨領域BNの3次元モデルである。図中の点SPは、顆間隆起の中央に位置し、人工関節置換手術を行なう場合にリーミング等を行なう基準点である。図29Bは、上記図29Aの内容に加え、抽出した2つに分かれた軟骨領域CT,CTの3次元モデルである。このように、ユーザが入力するのは、対象部位を囲む注目領域、Sagittal面とCoronal面に対する骨の概形の2項目とし、ユーザによる習熟度の影響を極力廃して操作入力の量を軽減させながら、安定して脛骨側の膝関節の軟骨部分の3次元画像を正確に抽出することが可能となる。処理時間1:骨領域の抽出処理時間と画素体積ボクセル数図30Aに骨領域の抽出処理時間と画素体積の測定結果を示す。体積は抽出処理で骨領域と判断された画素とした。テストデータ20例による軟骨領域の抽出処理時間は、画素体積ボクセル数に概ね比例する。テストに使用したPC10は、IntelCore2Duo登録商標、24GHzのCPU11と、2GBのメインメモリ13を搭載したものであり、大腿骨で1~2分、脛骨で3~9分程度の処理時間を要する。脛骨の処理時間が大腿骨の場合の3~5倍程度要しているが、これは動的輪郭法の処理が2方向分あり、1点あたりの処理データ量が15倍に増えることと、異方性拡散フィルタの3次元処理が2次元処理の時よりも1点当たりの処理データ量が3倍程度多いためである。2:軟骨領域の抽出処理時間と画素体積ボクセル数図30Bに軟骨領域の抽出処理時間と体積の測定結果を示す。体積は骨領域同様抽出処理で軟骨領域と判断された画素数とした。軟骨抽出は、大腿骨と脛骨の処理内容が異なるため、処理時間の長短を単純に比較することが難しい。表面積は、注目領域の対象となる面を和したものなので、実際に処理された表面積よりも大きく見積もられている。大腿骨は前後と下面の3面積の和、脛骨は上面の面積を採用している。このため、基本的に処理時間は体積に比例している。以上詳述した如く本実施形態によれば、ユーザによる習熟度の影響を極力廃して操作入力の量を軽減させながら、安定して関節の軟骨部分の3次元画像を正確に抽出することが可能となる。また特に上記実施形態では、Canny法で得られた画像上の輪郭に対して、動的輪郭法を適用し、膝関節骨領域の輪郭候補を得ることにしたので、比較的輝度変化の少ない輪郭部にも可変モデルの形状が適合できる。なお、上記実施形態は膝関節における大腿骨側と脛骨側それぞれの処理について説明したが、本発明は関節位置を特定するものではなく、その他の関節における軟骨領域の抽出にも適用可能であることは勿論である。その他、本発明は上述した実施形態に限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で種々に変形することが可能である。また、上述した実施形態で実行される機能は可能な限り適宜組み合わせて実施しても良い。上述した実施形態には種々の段階が含まれており、開示される複数の構成要件による適宜の組み合せにより種々の発明が抽出され得る。例えば、実施形態に示される全構成要件からいくつかの構成要件が削除されても、効果が得られるのであれば、この構成要件が削除された構成が発明として抽出され得る。以下に、本願出願の当初の特許請求の範囲に記載された発明を付記する。請求項1記載の発明は、関節領域のMRI3次元スキャン画像データを記憶した画像記憶部を備えたコンピュータが実行するプログラムであって、上記画像記憶部で記憶した画像データの関節骨領域を含む範囲を指定する範囲指定処理、上記指定した範囲内の関節骨領域の輪郭情報を抽出する輪郭抽出処理、同範囲内の関節骨領域を取得する領域拡張処理、及び、関節骨領域の輪郭情報に対する変形を受付ける可変モデル処理により骨領域を抽出する骨領域抽出ステップと、上記骨領域抽出ステップで抽出した骨領域から軟骨領域のエッジ位置群を検出するエッジ検出処理、及び、検出したエッジ位置群の連結から軟骨領域の境界を推定する最適経路検出処理により軟骨領域を抽出する軟骨領域抽出ステップとをコンピュータに実行させることを特徴とする。請求項2記載の発明は、上記請求項1記載の発明において、上記骨領域抽出ステップの可変モデル処理において、動的輪郭法により関節骨領域の輪郭情報を得ることを特徴とする。請求項3記載の発明は、上記請求項1記載の発明において、上記骨領域抽出ステップの輪郭抽出処理において、異方性拡散フィルタを用いて画像上のノイズを低減する画像平滑化処理を実行することを特徴とする。請求項4記載の発明は、上記請求項1記載の発明において、上記骨領域抽出ステップの領域拡張処理において、FastMarching法を用いることを特徴とする。請求項5記載の発明は、上記請求項1記載の発明において、上記骨領域抽出ステップの領域拡張処理において、動的輪郭法による輪郭曲線とFastMarching法による領域拡張で得られた領域との論理積で関節骨領域を取得することを特徴とする。請求項6記載の発明は、上記請求項1記載の発明において、上記軟骨領域抽出ステップの最適経路検出処理において、Dijkstraのアルゴリズム及びFastMarching法のすくなとも一方に基づいて軟骨の表面位置の候補となるエッジ点群から表面位置を推定することを特徴とする。請求項7記載の発明は、関節領域のMRI3次元スキャン画像データを記憶した画像記憶部と、上記画像記憶部で記憶した画像データの関節骨領域を含む範囲を指定する範囲指定処理、上記指定した範囲内の関節骨領域の輪郭情報を抽出する輪郭抽出処理、同範囲内の関節骨領域を取得する領域拡張処理、及び、関節骨領域の輪郭情報に対する変形を受付ける可変モデル処理により骨領域を抽出する骨領域抽出部と、上記骨領域抽出部で抽出した骨領域から軟骨領域のエッジ位置群を検出するエッジ検出処理、及び、検出したエッジ位置群の連結から軟骨領域の境界を推定する最適経路検出処理により軟骨領域を抽出する軟骨領域抽出部とを備えたことを特徴とする。"