JPWO2002087444A1 - ３次元関節構造計測方法 - Google Patents

Abstract

関節、関節軟骨の破壊部を再現性良く、正確かつ高速に抽出でき、関節及び破壊部を定量的かつ簡便に３次元解析できる方法であり、被検関節断面の２値化画像の関節内部空洞領域である骨髄腔領域を膨張収縮法を用いて穴埋め処理し、関節断面画像の外輪郭を抽出するステップにて生成される被検関節断面２値化画像、または前処理を行わない被検関節断面の２値化画像を積み上げた３次元画像にて３次元ラベリング処理を施し、評価対象となる関節画像を定義する３次元画像処理方法を提供する。

Description

技術分野
本発明は、関節断面撮影画像に基づいて関節の形態あるいは表面構造を測定する計測方法に関し、更に詳しくは、関節断面の２値化画像に基づいて関節から関節破壊部位（以下、破壊部）を抽出する計測方法において、関節から破壊部を自動的かつ高速かつ再現性よく抽出でき、非破壊的に３次元関節構造の解析の基礎となる関節と破壊部に分離した３次元画像が得られる計測方法に関する。更に上記計測方法によって得られた関節及び破壊部の３次元構造などの解析指標となる構造学的な種々のパラメータについて計測を実施できる方法に関する。
背景技術
従来、慢性関節リウマチや変形性関節炎等の関節炎における関節の破壊を評価する方法としては、習熟者による透過Ｘ線写真の読影を行なう方法が用いられている。透過Ｘ線写真の読影においては、ラーセン法（Ｌａｒｓｅｎ Ａ．ｅｔ ａｌ；Ａｃｔａ．Ｒａｄｉｏｌ．Ｄｉａｇ．１８．４８１（１９７７）、Ｌａｒｓｅｎ Ａ．ｅｔ ａｌ；Ｊ．Ｒｈｅｕｍａｔｏｌ．２２．１９７４（１９９５））や改良シャープス法（Ｓｈａｒｐ ＴＪ．ｅｔ ａｌ；Ａｒｔｈｒｉｔｉｓ Ｒｈｅｕｍ．２８．１３２６（１９８５）、ｖａｎ ｄｅｒ Ｈｅｉｊｄｅ ＤＭ．ｅｔ ａｌ；Ｂａｉｌｌｉｅｒｅｓ Ｃｌｉｎ．Ｒｈｅｕｍａｔｏｌ．１０．４３５（１９９６））、あるいは、ケルグナン・ローレンス法（Ｋｅｌｌｇｒｅｎ ＪＨ．ａｎｄ Ｌａｗｒｅｎｃｅ ＪＳ．ｅｔ ａｌ；Ａｎｎａｌｓ．Ｒｈｅｕｍａｔｉｃ．Ｄｉｓｅａｓｅｓ．１６．４９４（１９５７））等によりスコア化されて比較検討されている。しかし、これらのスコア化による評価方法は、基本的には定性的な評価であること、読影者の習熟度に大きく影響されること、再現性に乏しいこと等、種々の問題点を有している。
また、透過Ｘ線写真では直接観察することの不可能な関節軟骨に関して、例えばラーセン法やケルグナン・ローレンス法では関節構成骨間間隙、主に大腿骨と脛骨との間隔の変化を関節間隙の狭小化（Ｊｏｉｎｔ Ｓｐａｃｅ Ｎａｒｒｏｗｉｎｇ；ＪＳＮ）として間接的に評価することによって代替している。
しかし、透過Ｘ線写真は検査対象関節の特定方向からの通過像であり、本来３次元構造物である関節の構造的情報、詳しく言えば、関節表面、関節破壊部及び関節間隙（Ｊｏｉｎｔ Ｓｐａｃｅ；ＪＳ）等の体積、表面積、空間的情報等の３次元構造情報を著しく欠落させている。また、透過Ｘ線写真の他にＭＲＩや超音波を用いて直接関節軟骨の破壊を評価する方法も利用されているが、基本的には２次元的評価であり、透過Ｘ線写真による評価の問題点を本質的には解決していない。
先述のとおり関節は３次元構造物であり、慢性関節リウマチや変形性関節炎等の関節炎における関節破壊も３次元的に発生することから、関節破壊の評価も３次元にて解析することが望まれ、３次元解析が可能となれば、透過Ｘ線写真等による間接的な評価であっても、前述の２次元での解析における問題点は解消されるものである。
これまでに我々は、３次元骨構造を解析する方法を開発し、腰椎椎体部位等の３次元骨構造解析の有用性について報告している（日骨形態計誌．９．９７（１９９９）、日骨形態計誌．１０．５３（２０００））。さらに、この３次元骨構造解析方法を用いて、関節の３次元構造解析を実施し、その３次元構造の特徴を評価しうることを報告している（日骨代謝学会誌．１７．５７（１９９９））。
しかし、かかる方法においては、関節構成骨の内部構造、詳しく言えば関節構成骨内の骨梁構造に影響を受けてしまうことから、関節表面及び関節破壊部の３次元的特徴を正確に評価することは難しいと言える。
関節破壊部位の評価方法としては、関節軟骨もしくは関節軟骨下骨の形状を直接計測する方法と、関節軟骨もしくは関節軟骨下骨から破壊部を抽出し、その抽出された破壊部の形状を計測する方法が考えられる。しかるに、これらの概念に関して、求められる要望を充足しうる評価方法は存在しない。前述の３次元骨構造解析方法においても、腰椎椎体等の骨構造を詳細に検討するために、種々の画像処理方法を用いているが（例えば、ＷＯ００／７４５６７号明細書を参照）、関節の破壊部位を抽出する為の画像処理方法は開発し得ていない。従って関節から破壊部を再現性良く、正確かつ高速に抽出でき、関節及び破壊部を定量的かつ簡便に３次元解析できる方法の開発が要望されている。
同様に現時点における関節軟骨破壊の代替指標であるＪＳＮに関しても、関節構成骨間の間隙を３次元的に抽出する為の画像処理方法は開発し得ていない。従って関節からＪＳ領域を再現性良く、正確かつ高速に抽出でき、当該領域を定量的かつ簡便に３次元解析できる方法の開発が要望されている。これを解決することが本発明の目的である。
発明の開示
上記目的は以下の本発明にて達成される。すなわち本発明は、被検関節断面の２値化画像を積み上げた３次元画像上にて、３次元画像処理手段を用いることにより、評価対象関節の３次元画像のみを簡便にかつ正確に抽出することができる。さらに本発明は、関節表面の形状的特徴を抽出する為に、関節断面画像において関節骨骨髄腔を穴埋め処理することを特長とする。かかる処理によって得られた被検関節断面像の２値化画像をもとに関節の３次元再構築を行い、その再構築された３次元画像の解析を行った際に、従来の３次元構造解析では成し得なかった関節の内部構造、詳しく言えば関節内骨梁構造に影響されること無く、関節の表面構造を解析することができる。さらに、膨張収縮を組み合わせた方法或いは任意の直径の球体を凹凸に沿って走査する方法を適用することによって、関節表面の凹領域部、言い換えれば破壊部とそれ以外の部位とを構成した関節の３次元画像から３次元的に抽出する画像処理を特長とする。このように本発明は、関節の表面構造及び関節の破壊を定量的かつ客観的に３次元評価可能とするものである。
また本発明は、解析対象に関節を構成する２つの関節構成骨を指定し、３次元的な膨張収縮を組み合わせた方法或いは任意の直径の球体を表面の凹凸に沿って走査する方法を適用することによって、関節構成骨間のＪＳ領域、換言すれば関節軟骨部を３次元的に抽出する画像処理を特長とする。このように本発明は、関節構成骨間のＪＳ領域を定量的かつ客観的に３次元評価可能とするものである。
また関節表面の凹凸、言い換えれば破壊部を３次元的に抽出する方法として、膨張収縮を用いたアルゴリズムと球体による表面走査のアルゴリズムを用いることによって、被検査関節の形状、部位等、詳しく言えば、股関節、指間関節、踵、椎体等の変化に応じた最適な解析手法の提供を可能とするものである。
発明を実施するための最良の形態
以下、ラットの膝関節部位に適用した実施例に基づいて、実施の形態と共に順を追って説明する。
被検関節であるラットの膝関節部位の断面は、前述のＷＯ００／７４５６７号明細書に記載の方法と同様に、本例では撮像手段に膝関節部位の表面構造を計測するのに十分な焦点寸法及び解像度をもつ微小フォーカスＸ線断層撮影装置（以下、μＸ線ＣＴ）を用いて撮影した。なお、被検骨の撮影には、この他高解像度Ｘ線装置または核磁気共鳴画像診断装置（ＭＲＩ）またはフィルムスキャナー、顕微鏡からの出力画像等の２次元情報を発生する装置等が適用可能である。
μＸ線ＣＴを用いて撮像した場合、条件によっては関節軟骨の可視化は難しいが、関節軟骨下骨を評価することが可能である。この場合、計測対象は、関節軟骨下骨であり、その表面形状、凹凸の状態及び関節軟骨下骨の破壊部についての情報を得ることができる。また、ＭＲＩ等、軟骨を鮮明に撮像可能な装置を用いた場合は、計測対象は、関節軟骨であり、その表面構造、凹凸の状態及び軟骨破壊部についての情報を得ることができる。本発明においては、撮像手段に応じて計測対象は変化するものの、同一のアルゴリズムを用いて対象の表面構造、凹凸の状態及び破壊部について計測可能である。本発明では、関節の画像を提供する被検査対象が骨か軟骨かにかかわらず、関節表面の構造を計測できることから、被検査対象を「関節部」、抽出された破壊部領域を「破壊部」と表記し、一般的な関節表面及び関節破壊部の構造解析方法を提供する。
撮像手段で得られた該被検関節部位の画像信号を以下の画像処理装置で処理して、破壊部の計測を行う。本例の画像処理装置１０は、図１に示すように画像処理専用の画像処理ボード１１（本例ではシャープセミコンダクタ社製、ＧＰＢ−Ｋ）並びに画像記憶用に大容量ハードディスク１３を備え、画像表示用の外部モニター１２を搭載したパーソナルコンピュータを用いた。そして、図２に示すフローチャートの処理手順を収納し、撮像手段から得られた被検関節部位の画像信号を処理して所定の方法で破壊部領域を自動的に抽出する。
処理に先立ち、前述のＷＯ００／７４５６７号明細書に記載の方法と同様にして、まず、被検対象のラットの大腿骨関節部位を、μＸ線ＣＴ装置によりその軸方向に所定間隔で連続した断層像を撮影して、画像処理装置１０へ取り込む。なお、この画像は、関節表面形状を観察するのに十分な程度（１０〜３０μｍ）の空間解像度が得られるものがよい。本例において用いたμＸ線ＣＴ画像は、縦５１２画素×横５１２画素×高さ（軸方向）２１０スライス（最大５１２スライスまで可能）であり、断面の１画素あたりの寸法は３７．６μｍ、スライス間隔３０μｍ、各画素の輝度値ＣＴ（濃度値とも言う）は２^８階調で表現されている。この画像を、図２のステップ１の原画像入力により画像処理装置１０に読み込む。
図２のステップ２の２値化ステップを実行し関節軟骨下骨（以下、関節部）を定義する。その際、２値化の方法は、以下の判別分析法を用いる。
図３Ａはラット大腿骨関節付近断面の撮影原画像であり、関節部と背景部を明確化するため該撮影画像を２値化（物体無０か有１）することにより、図３Ｂを得る。実際の画像データは８ビットのデータであり、２値化することにより関節部データと背景データに２極化する必要がある。
画像の２個化方法の手段としては、３次元の画像データに対し判別分析法を使用した。判別分析法は、グループ内の濃度差が小さく、しかもグループ間の偏差が大きい場合に２値化の閾値を決定する際に使用することができる。式１〜式３を用いて、式１の分散比（Ｆ_０）が最大になるような閾値を求めることにより関節部を検出することができる。

以上の方法により該閾値を決定する方法を判別分析法という。この判別分析法にて検出した関節断面画像データをスライス方向に積み重ねれば関節部の３次元画像（図４）が得られる。
次に図２のステップ３の骨髄腔を抽出するステップを説明する。
３次元関節構造解析において実際に解析対象となるのは、骨表面の凹構造部、すなわち関節破壊部である。その場合、図３Ｂに示される複雑な関節内部構造が関節表面の構造的特徴を解析する上でノイズとなる為、あらかじめ関節内部空洞（以下、骨髄腔）を抽出しておき、関節内部構造の穴埋め画像、換言すれば関節部と骨髄腔部との和（ＯＲ）画像を作成する。この画像処理によって関節表面の構造的特徴を鋭敏に解析することが可能となる。
手順としては、前述の判別分析法にて２値化された関節部画像（図５Ｂ）に対してｎ回膨張処理（図５Ｃ）を施した後、最大穴部分を除いた穴部画像（図５Ｄ）を検出して、該穴部画像（図５Ｄ）と該膨張画像（図５Ｃ）の和（ＯＲ）画像（図５Ｅ）を得る。次にｎ回画像収縮処理（図５Ｆ）を施す。本実施例の場合には、該分断部分の間隙幅からｎ＝１０とした。
次に、該画像収縮画像（図５Ｆ）から該２値化画像（図５Ｂ）の差画像を求めれば骨髄腔画像（図５Ｇ）を検出することが可能である。ここで検出した骨髄腔画像データをスライス方向に積み重ねて、３次元上にてｍ回膨張処理、ｍ回収縮処理することにより骨髄腔の３次元画像（図６）が得られる。
次に図２のステップ４の３次元ラベリングを行うステップについて説明する。
該ラベリング処理は、関節部３次元画像データに対して、骨髄腔部を穴埋め処理した画像（図４と図６の和（ＯＲ）画像）を用いて実施し、図２のステップ５に示した解析対象関節部の抽出する為の前処理として行う。また、ここでいうラベリングとは、３次元画像データに対し指定された近傍同士に物体が存在する場合は、物体が連結されていると判断し、３次元上で連結されている物体にラベル番号を定義し、連結個数を数え上げるものである。
走査方法は図７左のように、要素（２，２，２）から要素（ｉ−１，ｊ−１，ｋ−１）まで走査する。その際、現在の要素をｘ０＝（ｉ，ｊ，ｋ）、その近傍で走査済みの要素（Ｘ１，Ｘ２・・・，Ｘ３）を図７右のように表し、さらに、Ｘｐのラベル（画像Ｌの値）をＬｐと省略する。連結が２６近傍の場合はＰ＝１３、連結が６近傍の場合はＰ＝３を用いる。
図８に、アルゴリズムのフローチャートを示す。ます、ステップ２にて、連結成分を表す変数λ（以下、λと称する）を０に初期化する。走査開始する要素を（２，２，２）とする。ステップ３にて現在の要素を（ｉ，ｊ，ｋ）と定義し、ステップ４にて入力画像（Ｆｉｊｋ）が０なのか１なのかで処理を分岐する。０の場合は、そこに物体が存在しないことになるのでラベル画像Ｌｉｊｋを０とし、ステップ８を実行する。一方、ステップ４にて入力画像（Ｆｉｊｋ）が１の場合は、ステップ５において｛Ｔ（Ｌｐ）：Ｐ＝０，１，２、・・・，１３｝の中に異なる正の値がｎ種類あるとし、それらを小さい順に、Ｌ１，Ｌ２，・・・，Ｌｎとする。
そしてステップ６にて、そのｎの値に対応してステップ７にて処理を分岐する。ｎ＝０の場合は、初めて出現する物体であるため、λ←λ＋１とし、ラベル表Ｔ（λ）←λを代入、ラベル表Ｌｉｊｋ←λを代入し、ステップ８を実行する。ｎ＝１の場合は、連結成分が同じになる物体が１つしかないことになるので、ラベル表Ｌｉｊｋ←λを代入し、ステップ８を実行する。それ以外の場合は、連結成分が同じになる物体が複数個存在することになるので、ラベル番号が一番小さいＬ１をラベル画像Ｌｉｊｋへ代入し、Ｔ（γ）＝Ｌｐ（２≦ｐ≦Ｐ、１≦γ≦λ）となるすべてのＴ（γ）に対してＴ（γ）＝Ｌ１を代入し、ステップ８を実行する。ステップ８では、走査する要素が残ってなければステップ９へ、残っていればステップ３へ戻り走査する要素がなくなるまでステップ３からステップ８を繰り返す。ステップ９にてラベル表Ｔ（λ）は歯抜けの番号が存在することになるので、通し番号になるよう割り当てなおす。ステップ１０にて、ラベル表Ｔ（λ）の格納されているラベル番号の最大値を検出し、それを最大ラベル番号：Ｌｍａｘとする。次にステップ１３にて、すべてのラベル画像を走査し、０より大きいならばステップ１２にて、通し番号にしたラベル表Ｔ（Ｌｉｊｋ）をＬｉｊｋに代入すれば、すべての処理が完了する。最終的なラベル画像：Ｌｉｊｋが３次元ラベリングデータを格納しており、Ｌｍａｘにその最大ラベル番号が格納されている。実際のラット膝関節部位を用いて２６近傍にてラベリングした結果を図９に示す。
次に図２のステップ５に示した解析対象関節部の抽出について説明する。該抽出処理は、通常、図９に示したように被検関節領域には、関節構成成分が多数存在する為、これらから解析対象を抽出する目的で行う。
図２のステップ４に示した３次元ラベリング処理によって、被検関節領域に存在する独立した全ての関節構成成分は、独立したラベル番号が与えられていることから、任意のラベル番号の３次元構造体、換言すれば、任意の関節構成成分を解析対象関節部に指定することが可能である。
また、本発明においては、膝関節部位を解析対象とした場合、各関節構成骨の重心位置の相対関係から評価対象となる関節画像すなわち大腿骨、脛骨、膝蓋骨、半月板を自動的に判別・同定し、抽出するアルゴリズムを開発した。
本実施例では、ラット膝関節部位（図９）を用いて解析対象となる大腿骨部、脛骨部、膝蓋骨部、半月板部を抽出する方法を説明する。抽出方法を説明するにあたり、撮影されるＣＴ断層スライス画像は脛骨から大腿骨方向に撮影された順でスライス画像群が構成されている（スライス画像番号が若い方が脛骨側）ものとする。まず脛骨部の抽出方法は、ＣＴ断層スライス画像の一番目に接していてボクセル個数の一番大きい特徴を利用して抽出する。ボクセル個数とは、ラベル画像Ｌｉｊｋにて各ラベル番号が存在する頻度をカウントすることにより各ラベル画像が持っているボクセルの個数として定義される。大腿骨部の抽出方法は、ＣＴ断層スライス画像の最終番目に接していてボクセル数の一番大きい特徴を利用して抽出する。膝蓋骨の抽出は、同じくＣＴ断層スライス画像の最終番目に接していてボクセル数の２番目に大きい特徴を利用して抽出する。半月板の抽出方法は、図１０に示すように各解析対象関節部の重心座標を計算し、膝蓋骨の重心座標と脛骨の重心座標から定義される領域の上部半分（ＣＴ断層スライス画像番号の若い方向）に存在するラベル群を半月板として定義することが可能である。抽出された大腿骨を値３２、脛骨を値６４、膝蓋骨を値９６、半月板を値１２８として定義することによって、各解析対象骨毎に以下に示す３次元構造解析を実施することが可能である。
次に図２のステップ６にて示した破壊部抽出処理について、ステップの手順を追って説明する。本発明において、凹領域の抽出に膨張収縮を利用する方法（以下、膨張収縮法）と、物体境界画素に対し球体を逐次当てはめていき、物体と球体との重なる割合により、その着目画素を抽出するかどうかを判定する方法（以下、球体当てはめ法）を開発した。まずは前者の膨張収縮法について説明する。図１１に簡単に原理を２次元にて示す。まずステップ１にて説明用の入力画像を用意した。この画像は物体が白色（以下、１−画素）、背景が黒色（以下、０−画素）とし、凹構造として２ｎピクセルの長方形の穴を用意した。ピクセルとは画像の大きさを表す単位である。ステップ２にて画像をｎ回膨張させた結果を示す。この場合、画像境界も大きくなるが２ｎピクセルの穴も埋まることになる。ステップ３にて画像をｎ回収縮させた結果を示す。画像全体の大きさは入力画像と一致するが凹領域は穴埋めされているのが分かる。ステップ４にて収縮画像（ステップ３）から入力画像（ステップ１）を差し引くと凹領域のみの画像が抽出可能になる。この膨張収縮法の特徴としては、利点として指定した膨張回数の２倍の凹領域が抽出可能で、比較的計算時間も早いことが上げられる。この膨張収縮法を用いて、関節破壊を有するラットの大腿骨関節部の破壊部を抽出した実施例を図１２に示す。
次に後者の方法である球体当てはめ法について手順を追って説明する。図１３に簡単に原理を示す。まずステップ１にて説明用の入力画像を用意した。この画像は物体が白色（以下、１−画素）、背景が黒色（以下、０−画素）とし、凹構造として２ｎピクセルの長方形の穴を用意した。次にステップ２にて物体の境界画素を抽出する。境界画素の抽出方法は、着目した１−画素において隣り合う近傍画素に０−画素が存在すれば、その着目画素は境界となる。次にステップ３で、着目する境界画素上に半径ｒ１の判定領域定義用球体（以下、球体Ｃ１）を当てはめ物体と重なり合う体積を数え上げる。この例の場合、領域（Ａ＋Ｂ＋Ｃ）がそれに相当する。また、ある程度局率を持たせる形で穴埋めを進行させたい場合は、半径ｒ２の局率定義用球体（以下、球体Ｃ２）を球体Ｃ１上に中心を置くように当てはめ球体Ｃ１と球体Ｃ２が重ならない領域（Ａ）を判定体積とする。この着目した１−画素を穴埋めするか、しないかは以下の式４にて判定する。
式４： 領域（Ａ＋Ｂ＋Ｃ）の体積 ＞ 領域（Ａ）の体積
領域（Ａ＋Ｂ＋Ｃ）の体積が、領域（Ａ）の体積よりも大きい場合、着目した境界画素である１−画素は、穴埋めを実施する処理とする。この処理を画像全体に対して穴埋めする場所がなくなるまで走査していけば、ステップ４の穴埋め画像（図１３Ｄにおける１３１）を定義可能となる。この球体当てはめ法の特徴としては、評価対象となる３次元構造体が有する表面の緩やかな凹凸構造に対して寛容であり、急激な表面構造の変化を選択的に抽出できることである。図７左図の走査の方法で、この球体当てはめ法を用いて関節破壊を有するラットの大腿骨関節部位の破壊部を抽出した実施例を図１４に示す。
次に図２のステップ７の破壊部構造パラメータについて計測方法を説明する。破壊部構造パラメータとして定義可能な項目を、図１５へ示す。
前提条件として関節部領域と骨髄腔領域と破壊部領域と背景領域のボクセル値の取り扱いを以下のように定義する。
（ａ）関節部領域 →１
（ｂ）骨髄腔領域 →２
（ｃ）破壊部領域 →３
（ｄ）背景領域 →０
ボクセルとは、図１６の升目の最小単位であり、ボクセル値とは、その升目上に定義されている値のことである。（ａ）は図２のステップ２にて２値化した関節部領域データに相当し、（ｂ）は図２のステップ３にて抽出した骨髄腔領域に相当し、（ｃ）は図２のステップ６にて抽出した破壊部領域に相当し、関節部領域のボクセル値を１、骨髄腔領域のボクセル値を２、破壊部領域のボクセル値を３、背景領域のボクセル値を０と定義する。
まず、図１５のパラメータ１、パラメータ２の関節体積、関節表面積について説明する。関節体積、関節表面積を計測するにあたっては、物体の境界面を定義する方法（改良マーチングキューブス法：日本特許登録番号０２７９８５６１号）を適用し、境界表面に三角形を連続する様に定義する画像処理手段を用いる。マーチングキューブス法とは、ボクセルを頂点とした２×２×２の８個のキューブのパターン（図１７）から考えられる２１通りのボクセル値のパターン上にて予め境界面を図１８のように定義し、実モデルとのパターンマッチングにて境界面を定義していく方法である。関節の境界を定義するには、ボクセル値が１の関節部領域とボクセル値が２の骨髄腔部を１つの物体と考え、それ以外の領域（ボクセル値が０または３）で二極化されるデータの境界に面を定義することである。そこで、関節部データより２×２×２のボクセル領域でサンプリングされたキューブと、図１８のキューブ上のパターンマッチングにより境界面を定義し、これらのキューブデータおよび境界面データを数え上げていくことにより、関節体積と関節表面積が計測可能になる。
次に図１５のパラメータ３である関節ＢＳ／ＢＶの求め方について説明する。関節ＢＳ／ＢＶについては、式５にて算出する。
式５： 関節ＢＳ／ＢＶ＝関節表面積／関節体積
関節ＢＳ／ＢＶは、単位体積あたりの表面積の量を表し、関節部表面構造の複雑性を表す指標の１つである。
次に図１５のパラメータ４である関節表面凹凸度（Ｓｕｒｆａｃｅ Ｉｒｒｅｇｕｌａｒｉｔｙ）の求め方について説明する。関節表面凹凸度については式６にて算出する。

図１９に示すように体積（Ｖ１）と表面積（Ｓ１）は、膨張前の体積、表面積を示し、体積（Ｖ２）と表面積（Ｓ２）は画像膨張後（１ボクセル分の膨張とする。）の体積、表面積を示す。この数式によって求められた値は、凸構造部位では正となり、凹構造部位では負となる。従って全体の構造として、凸構造が多数を占める構造体では正の値となり、凹構造が多数を占める構造体では負の値となり、凸構造と凹構造がほぼ等しい場合、ゼロ近傍の値となる。
次に図１５のパラメータ５、パラメータ６の破壊部体積、破壊部表面積について説明する。破壊部体積、破壊部表面積を計測するにあたっては、図１５のパラメータ１およびパラメータ２と同様に物体の境界面を定義する方法（改良マーチングキューブス法：日本特許登録番号０２７９８５６１号）を適用し、境界表面に三角形を連続する様に定義する画像処理手段を用いることにより破壊部体積、破壊部表面積が求まる。
次に図１５のパラメータ７である破壊部内側表面積の求め方について説明する。破壊部内側表面積とは関節部領域（ａ）と、破壊部領域（ｃ）が重なっている部分の表面積を示す。この場合、関節部領域（ａ）と破壊部領域（ｃ）のみがキューブ上に存在するパターンのみに対しマーチングキューブス法を適応することにより、破壊部内側表面積が求められる。破壊部内側表面積は換言すれば、関節表面における破壊部の表面積を意味しており、関節破壊の程度を判断する上で重要な指標の１つである。
次に図１５のパラメー８である破壊部個数の求め方について説明する。破壊部領域（ｃ）に対し、図２のステップ４で説明した３次元ラベリングを適応することで、画像上で独立している固体の個数が求められ、この個数が破壊部個数に相当する。
次に図１５のパラメータ９である破壊部ＢＳ／ＢＶについては、式７にて求められる。

破壊部ＢＳ／ＢＶは、単位体積あたりの表面積の量を表し、破壊部構造の複雑性を表す指標の１つである。
次に図１５のパラメータ１０である破壊部体積率の求め方について説明する。破壊部体積率の定義としては、式８を計算することによって求められる。
式８： 破壊部体積率＝破壊部体積／（骨体積＋骨髄腔体積＋破壊部体積）
次に図１５のパラメータ１１である破壊部内側表面積率の計算方法について説明する。破壊部内側表面積率の定義としては、式９を計算することによって求められる。
式９： 破壊部内側表面積率＝破壊部内側表面積／｛（骨領域ＯＲ骨髄腔領域）の表面積｝
次に図１５のパラメータ１２である破壊部平均体積の計算方法について説明する。破壊部平均体積とは、空間上で独立している破壊部の１個あたりの平均体積をあらわす。式１０にその計算方法を示す。
式１０： 破壊部平均体積＝破壊部体積／破壊部個数
次に図１５のパラメータ１３である破壊部平均表面積の計算方法について説明する。破壊部平均表面積とは、空間上で独立している破壊部の１個あたりの平均表面積をあらわす。式１１にその計算方法を示す。
式１１： 破壊部平均表面積＝破壊部表面積／破壊部個数
数７〜数１１で求められるパラメータ９〜パラメータ１３までの各構造パラメータは、関節破壊部の３次元構造体としての構造的特徴を示すものであり、関節破壊部の状態を判断する上で重要な指標群である。
次に図２のステップ８にて示したＪＳ領域抽出について、ステップの手順を追って説明する。
ＪＳ領域は対になる関節構成骨の間隙であり、関節破壊部領域と同様にＪＳ領域に関しても解析対象となる関節構成骨に対して膨張収縮法、球体当てはめ法を適応することによりＪＳ領域の３次元画像データの抽出が可能である。
前述の図２のステップ４に示した３次元ラベリング処理によって、本発明においては、被検関節領域に存在する独立した全ての関節構成成分は、独立したラベル番号が与えられている。これにより任意のラベル番号を有する２つの３次元構造体、換言すれば、任意の対となる関節構成骨を解析対象に指定することが可能である。また、同様に本発明においては、膝関節部位を解析対象とした場合、各関節構成骨の重心位置の相対関係から評価対象となる関節画像すなわち大腿骨、脛骨、膝蓋骨、半月板を自動的に判別し、ラベル番号を与えることができる。これにより、膝の対となる関節構成骨、つまり大腿骨・脛骨を解析対象として自動的に指定することが可能である。
ＪＳ領域の定義方法のフローチャートを図２０に示し、ラット膝関節つまり大腿骨と脛骨のモデルを実施例として具体的なステップを説明する。
まず図２０のステップ１にて解析対象となる２つの関節構成骨の合成画像（画像１）を作成する。次に図２０のステップ２にて画像１に対し関節破壊部領域と同様に膨張収縮法（または球体当てはめ法）を実施し、画像１との差画像（画像２：図２１Ａの２１３）を作成する。次に図２０のステップ３にて大腿骨のみの画像で膨張収縮法を実施し大腿骨の関節破壊部領域を抽出して画像３（図２１Ｂの２１４）を作成する。次に図２０のステップ４にて脛骨のみの画像で膨張収縮法を実施し脛骨の関節破壊部領域を抽出して画像４（図２１Ｃの２１５）を作成する。次に図２０のステップ５にて画像２から画像３と画像４の差画像を作成すると、ＪＳ領域（画像５：図２１Ｄの２１６）を定義することが可能である。ここで示した膨張収縮法を球体当てはめ法に置き換えてもＪＳ領域の抽出は可能である。
次に図２のステップ９のＪＳ領域パラメータについて計測方法を説明する。ＪＳ領域パラメータとして定義可能な項目を図２２へ示す。前提条件として関節部領域と骨髄腔領域とＪＳ領域と背景領域のボクセル値の取り扱いを以下のように定義する。
（ｅ）大腿骨関節部領域 →１
（ｆ）脛骨関節部領域 →２
（ｇ）骨髄腔領域（大腿骨側）→３
（ｈ）骨髄腔領域（脛骨側） →４
（ｉ）ＪＳ領域 →５
（ｊ）背景領域 →０
まず、図２２のＪｏｉｎｔ Ｓｐａｃｅ Ｖｏｌｕｍｅ（以下ＪＳＶと呼ぶ）について説明する。ＪＳＶとはＪＳ領域の体積のことである。ＪＳＶを計測するにあたっては、マーチングキューブ法を適用する。具体的にはＪＳ領域（ｉ）とそれ以外の領域（ｅ）、（ｆ）、（ｇ）、（ｈ）、（ｊ）に対して、２×２×２のボクセル領域でサンプリングされたキューブと、図１８のキューブ上のパターンマッチングにより境界面を定義し、これらのキューブデータおよび境界面データを数え上げていくことにより、ＪＳＶが計測可能になる。
次に図２２のＪｏｉｎｔ Ｓｐａｃｅ Ｍｉｎｉｍｕｍ Ｄｉｓｔａｎｃｅ（以下ＪＳＭＤと呼ぶ）の計算方法について説明する。ＪＳＭＤとは、対となる関節構成骨間の距離の最短長、この例では大腿骨関節部と脛骨関節部間の最短距離のことである。ＪＳＭＤの計算方法のフローチャートを図２３に示した。
図２３のステップ１にてＪＳ領域（ｉ）と大腿骨関節部領域（ｅ）および大腿骨側骨髄腔領域（ｇ）が接する大腿骨側のボクセル座標をリスト化（リスト１）する。次に図２３のステップ２にてＪＳ領域（ｉ）と脛骨関節部領域（ｆ）および脛骨側骨髄腔領域（ｈ）が接する脛骨側のボクセル座標をリスト化（リスト１）する。以上の処理で、図２７に示すようにＪＳ領域に接する大腿骨側のボクセルリスト（２４１）および脛骨側のボクセルリスト（２４２）が抽出できたことになる。図２３のステップ３にて大腿骨側のボクセル座標リスト（２４１）から脛骨側のボクセル座標リスト（２４２）へ距離データを求め最短の距離をＪＳＭＤとする。
次に本発明において解析可能となった３次元構造パラメータについて解析実施例を示し、その意義・有用性について説明する。
図２５に正常ラットの大腿骨関節部（以下、正常）と関節破壊を有するラットの大腿骨関節部（以下、病態）について、図１５で示した破壊部に関するパラメータについて解析した結果を示す。
関節部（この実施例では関節軟骨下骨）を直接計測したパラメータ（図２５Ａ〜図２５Ｄ）においても、抽出した破壊部を計測したパラメータ（図２５Ｅ〜図２５Ｍ）においても、正常と病態の間で破壊部の増加を検出することができた。
病態における関節ＢＳ／ＢＶの増加（図２５Ｃ）は、表面構造の複雑性の増加を意味し、関節表面凹凸度の低下（図２５Ｄ）は、表面凹構造の増加を意味し、どちらのパラメータの変化も病態において関節の破壊が進行したことを意味している。特に関節表面凹凸度に関しては、本特許において実施した関節内部の穴埋め処理をすることによってはじめて関節表面の形状変化を特異的に抽出できたものである。比較として関節内部の穴埋め処理を実施せずに関節表面凹凸度を計測した結果を図２６に示した。穴埋め処理を実施した場合（図２５Ｄ）と異なり、病態に対して正常が低値となり、正常の方が病態に対して表面凹構造が多いこと示している。これは病態対して正常の方が、関節内部構造体すなわち骨梁が多く、複雑な構造をしていることに起因する。これらのことから、本発明によって関節表面の構造的変化の検出感度が飛躍的に向上したと言える。
さらに２種類のアルゴリズム（膨張収縮法、球体当てはめ法）によって抽出された破壊部の計測結果（図２５Ｅ〜図２５Ｍ）においては、いずれのパラメータにおいても正常に対して病態で増加しており、本発明により病態において関節の破壊が増加していることを検出できた。また、膨張収縮法に対して球体当てはめ法は、低値をとることが見てとれるが（図２５Ｅ〜図２５Ｍ）、これは上記２種類の抽出方法の特徴に起因する。更に詳しく言えば、膨張収縮法は表面構造の微細な変化まで鋭敏に抽出するのに対して、球体当てはめ法は、表面の緩やかな凹凸に対して寛容である。目的とする評価対象に応じて２種類のアルゴリズム（膨張収縮法、球体当てはめ法）を区別して用いることによって、より正確かつ詳細な表面構造の変化を検出できるものである。
次に正常モルモットの膝関節部（以下、正常）と関節破壊を有するモルモットの膝関節部（以下、病態）において、ＪＳ領域に関して解析した結果について説明する。
図２７に正常と病態に関して、膨張収縮法によって抽出したＪＳ領域を示した。両例ともにＪＳ領域のみが正確に抽出されることが確認され、本発明によって、これまで困難であったＪＳ領域の３次元的な認識及び抽出が可能であることが示された。当該ＪＳ領域は球体当てはめ法によっても同様に抽出可能である。
次に図２８に図２７で示した正常及び病態のＪＳ領域に対してＪＳＭＤを計測した結果を示した。病態においてＪＳＭＤの減少が認められ、正常と比較して関節軟骨の破壊もしくは消失が生じていることが観察された。
以上のことから、本発明によってＪＳ領域の３次元的な視覚認識と３次元的な構造パラメータを客観的かつ定量的に再現性良く解析することが可能となった。
発明の効果
本発明では、関節内部構造の穴埋め処理をすることによって関節表面の構造を感度良く抽出でき、さらに２種類の画像処理方法の開発により、関節部から凹構造領域すなわち破壊部を高速かつ簡単に再現性良く抽出することが可能となった。さらに得られた関節表面および破壊部の画像データは、ともに３次元構造解析を実施することによって、関節破壊に関する構造的な指標を提供することができる。また、関節を構成する対となった関節構成骨を解析対象とすることによって、ＪＳ領域の３次元構造を高速かつ再現性良く解析することが可能となった。
また、同一のアルゴリズムを用いて関節表面及び破壊部の抽出をすることが可能であることから、複数の被解析対象関節を解析する場合、すなわち同一患者の経時的変化の解析や薬剤の投与・非投与による治療効果の判定等においても、観察者の主観が介入することなく関節表面及び破壊部の構造的特性を再現性良く、客観的かつ定量的に評価することができる。
さらに、本発明において定義された関節部および関節破壊部を可視化ツールを用いて表示することが可能であり、関節表面構造および関節部に対する関節破壊部の３次元空間的情報を視覚的に観察者に提供できるものである。また、本発明による３次元関節構造解析は、膝関節のみでなく他の関節部位への適用についても可能である。
さらに本発明は、μＸ線ＣＴ撮影装置をはじめ臨床用ＣＴ撮影装置、ｐＱＣＴ撮影装置またはＭＲＩ等から得られる診断画像情報全般を関節の３次元構造計測に応用できると期待できる。
【図面の簡単な説明】
図１は、本発明を実施する装置例の構成の説明図である。
図２は、関節部、関節破壊部及びＪＳ領域を計算するための全体フローチャートである。
図３は、ラット大腿骨関節付近断面の原画像（図３Ａ）と２値化画像（図３Ｂ）である。
図４は、ラット大腿骨関節軟骨下骨の３次元可視化画像である。
図５（図５Ａ〜図５Ｇ）は、ラット大腿骨関節断面から骨髄腔を抽出するステップを示す画像群である。
図６は、ラット大腿骨関節付近の骨髄腔の３次元可視化画像である。
図７は、３次元ラベル付け手順の説明図である。
図８は、３次元ラベル付け手順のフローチャートである。
図９は、ラット膝関節部位の関節構成骨の３次元表示である。
図１０は、ラット膝関節部位における各関節構成骨の重心位置の３次元表示である。
図１１（図１１Ａ〜図１１Ｄ）は、膨張収縮法を用いた破壊部領域の抽出方法の手順である。
図１２は、膨張収縮法を用いて抽出した破壊部領域の３次元表示である。
図１３（図１３Ａ〜図１３Ｄ）は、球体当てはめ法を用いた破壊部領域の抽出方法の手順である。
図１４は、球体当てはめ法を用いて抽出した破壊部領域の３次元表示である。
図１５は、破壊部計測で使用するパラメータリストである。
図１６は、ボクセル定義を説明する図である。
図１７は、マーチングキューブス法のキューブ定義を説明する図である。
図１８は、マーチングキューブス法で使用するキューブパターンと境界を示す図である。
図１９は、骨表面凹凸度（Ｓｕｒｆａｃｅ Ｉｒｒｅｇｕｌａｒｉｔｙ）の計算方法を示す図である。
図２０は、ＪＳ領域の抽出方法の手順である。
図２１（図２１Ａ〜図２１Ｄ）は、ＪＳ領域を定義するための説明図である。
図２２は、ＪＳ領域パラメータで使用するパラメータリストである。
図２３は、Ｊｏｉｎｔ Ｓｐａｃｅ Ｍｉｎｉｍｕｍ Ｄｉｓｔａｎｃｅを計算する手順である。
図２４は、Ｊｏｉｎｔ Ｓｐａｃｅ Ｍｉｎｉｍｕｍ Ｄｉｓｔａｎｃｅを計算する手順の説明図である。
図２５は、図１５で示したパラメータに関するラット大腿骨関節軟骨下骨の計測結果（図２５Ａ〜図２５Ｍ）である。
図２６は、関節内部構造の穴埋め処理を実施しなかった場合の表面凹凸度の計測結果である。
図２７は、膨張収縮法を用いて抽出したモルモットの膝関節のＪＳ領域（図２７Ａ、図２７Ｂ）である。
図２８は、図２７で示したモルモットの膝関節ＪＳ領域のＪＳＭＤの計測結果である。
図中の符号は、それぞれ以下のことを表わす。
１０ 画像処理装置
１１ パーソナルコンピュータ
１２ 画像処理ボード
１３ ハードディスク
９１ 脛骨
９２ 大腿骨
９３ 膝蓋骨
９４ 半月板
１０１ 脛骨の重心
１０２ 膝蓋骨の重心
１０３ 半月板の重心
１３１ 球体当てはめ法によって抽出された凹部領域
２１１ 脛骨
２１２ 大腿骨
２１３ 膨張収縮法にて抽出されたＪＳ領域及び脛骨・大腿骨の関節破壊部位
２１４ 大腿骨側の関節破壊部位
２１５ 脛骨側の関節破壊部位
２１６ ＪＳ領域
２４１ ＪＳ領域に接する大腿骨側のボクセルリスト
２４２ ＪＳ領域に接する脛骨側のボクセルリスト

Claims (22)

1. 被検関節断面の２値化画像の関節内部空洞領域である骨髄腔領域を膨張収縮法を用いて穴埋め処理し、関節断面画像の外輪郭を抽出するステップを有する画像処理方法。
2. 請求項１の画像処理方法にて生成される被検関節断面の２値化画像または前処理を行わない被検関節断面の２値化画像を積み上げた３次元画像にて３次元ラベリング処理を施し、評価対象となる関節画像を定義する３次元画像処理方法。
3. 請求項２の画像処理方法にて作成される関節部位の３次元ラベル画像より、各関節構成骨の重心位置の相対関係から評価対象となる関節画像を判別・同定し、抽出する画像処理方法。
4. 請求項２の画像処理方法にて作成される膝関節位の３次元ラベル画像より、各関節構成骨の重心位置の相対関係から評価対象となる大腿骨、脛骨、膝蓋骨、及び／または半月板の関節画像を判別・同定し、抽出する画像処理方法。
5. 請求項２の画像処理方法にて生成される関節骨３次元画像に対し、３次元画像膨張法と３次元画像収縮法を用いて、関節破壊部位の３次元画像を抽出するステップを有する画像処理方法。
6. 請求項２の画像処理方法にて生成される関節３次元骨画像に対し、３次元球体当てはめ法を用いて、関節破壊部位の３次元画像を抽出するステップを有する画像処理方法。
7. 請求項３の画像処理方法にて生成される関節骨３次元画像に対し、３次元画像膨張法と３次元画像収縮法を用いて、関節破壊部位の３次元画像を抽出するステップを有する画像処理方法。
8. 請求項３の画像処理方法にて生成される関節３次元骨画像に対し、３次元球体当てはめ法を用いて、関節破壊部位の３次元画像を抽出するステップを有する画像処理方法。
9. 請求項４の画像処理方法にて生成される関節骨３次元画像に対し、３次元画像膨張法と３次元画像収縮法を用いて、関節破壊部位の３次元画像を抽出するステップを有する画像処理方法。
10. 請求項４の画像処理方法にて生成される関節３次元骨画像に対し、３次元球体当てはめ法を用いて、関節破壊部位の３次元画像を抽出するステップを有する画像処理方法。
11. 請求項２の画像処理方法にて生成される関節骨３次元画像について、表面凹凸度（Ｓｕｒｆａｃｅ Ｉｒｒｅｇｕｌａｒｉｔｙ）を計測することを特徴とする関節構造評価方法。
12. 請求項１〜１０の画像処理方法にて生成される３次元構造体について、体積（ＢＶ）、表面積（ＢＳ）、ＢＳ／ＢＶ、個数、平均体積（体積／個数）、及び／または平均表面積（表面積／個数）の構造パラメータを計測することを特徴とする関節構造評価方法。
13. 請求項２の画像処理方法にて生成される関節３次元骨画像の内、任意の２つの解析対象に対して、３次元画像膨張法と３次元画像収縮法を用いて、３次元間隙画像を抽出するステップを有する画像処理方法。
14. 請求項２の画像処理方法にて生成される関節３次元骨画像の内、任意の２つの解析対象に対し、３次元球体当てはめ法を用いて、３次元間隙画像を抽出するステップを有する画像処理方法。
15. 請求項３の画像処理方法にて定義される関節３次元骨画像の内、任意の２つの解析対象に対して、３次元画像膨張法と３次元画像収縮法を用いて、３次元関節間隙画像を抽出するステップを有する画像処理方法。
16. 請求項３の画像処理方法にて定義される関節３次元骨画像の内、任意の２つの解析対象に対して、３次元球体当てはめ法を用いて、３次元関節間隙画像を抽出するステップを有する画像処理方法。
17. 請求項４の画像処理方法にて定義される関節３次元骨画像の内、任意の２つの解析対象に対して、３次元画像膨張法と３次元画像収縮法を用いて、３次元関節間隙画像を抽出するステップを有する画像処理方法。
18. 請求項４の画像処理方法にて定義される関節３次元骨画像の内、任意の２つの解析対象に対して、３次元球体当てはめ法を用いて、３次元関節間隙画像を抽出するステップを有する画像処理方法。
19. 請求項１３〜１８の画像処理方法にて生成される３次元構造体について、Ｊｏｉｎｔ Ｓｐａｃｅ Ｖｏｌｕｍｅ（ＪＳＶ）、Ｊｏｉｎｔ Ｓｐａｃｅ Ｍｉｎｉｍｕｍ Ｄｉｓｔａｎｃｅ（ＪＳＭＤ）の構造パラメータを計測することを特徴とする関節構造評価方法。
20. 請求項１〜１９記載のいずれかの方法を用いて、関節関連疾患動物モデルにおいて関節の構造を評価する評価方法。
21. 請求項１〜１９記載のいずれかの方法を用いて、関節関連疾患患者の関節の構造を評価する評価方法。
22. 請求項１〜１９記載のいずれかの方法を用いて、関節関連疾患患者の関節破壊を診断する診断方法。

Images