JP2010527647A - 半自動式輪郭検出方法 - Google Patents

Abstract

構造物を含む画像を処理することによって、その画像内の構造物の輪郭の位置を特定する方法が提供される。画像内の構造物は、その上に３乃至１０個のランドマークポジションがアノテートされており、その構造物を含む画像を表現するデジタルデータの開始集合を用意するステップが実行される。そして、構造物の統計モデルを、画像上にアノテートされたランドマークポジションにフィットさせ、構造物の輪郭の初期推定値を生成するステップが実行される。更に、推定された輪郭に隣接するポイントから得られた階調レベル情報を使用して、反復的に輪郭境界線を移動させて構造物の輪郭の最終推定値を生成するステップが実行される。

Description

本発明は、画像内の構造物の輪郭を位置特定する方法に関する。

骨粗鬆症（osteoporosis）は、骨ミネラル密度（bone mineral density：ＢＭＤ）が減少し、骨の微細構造が破壊され、かつ骨に含まれる非コラーゲン性タンパク質の量と種類が変化する骨疾患である。この骨疾患に冒された骨はより骨折しやすくなっている。骨粗鬆症は、世界保健機関（World Health Organisation：ＷＨＯ）によれば、二重エネルギーＸ線吸収測定法（dual X-ray absorptiometry：ＤＸＡ）で測ったときの骨ミネラル密度が（２０歳の性別が同じ健常者の平均の）最大骨量から２．５ＳＤ（標準偏差）以上減少した状態か、あるいはその他の骨折しやすい状態として定義されている。そのホルモン成分のせいで、女性、特に閉経後の女性は、男性よりもこの疾患にかかりやすい。

骨粗鬆症性の骨折は、骨粗鬆症に罹っていない人なら通常骨折しないような軽度のストレスのもとで起こる骨折である。一般的な骨折は、脊柱、臀部および手首で起こる。椎骨骨折は骨粗鬆症でもっともよく起こる骨折である。若年患者で起これば、それらは脊柱および臀部の将来起こりえる骨折のリスクの良い指標である。これら２つはもっとも重篤なケースで、移動が制限されるか、場合によっては身体障害になる。特に臀部骨折は、通常、すぐに手術が必要であるが、これには、例えば、深部静脈血栓症または肺塞栓症などの他の重大なリスクが伴う。骨粗鬆症患者は骨折の合併症により死亡率が高まるが、ほとんどの患者は骨折で死ぬよりもそのような二次的な病気で死亡する。

椎骨骨折は、従来的にラテラルＸ線により検出され、等級分けされる。放射線科医師による主観的な画像診断とは別に、標準的な６点形態計測（six-point morphometry）が一般的に利用される。この技術では、６個のランドマーク（目印）が両方の椎骨端板（vertebra endplate）の角部（corners）と中間点に配置され、前高（anterior height）、中高（middle height）、後高（posterior height）が指定される。これらの測定値は骨折等級を計算するために使用することができる。現行の臨床試験では、骨折した椎骨は３つの高さの１つが他のどの高さよりも少なくとも２０％大きいものとして定義される。

６点表現は、画像内の重要な情報のほとんどを捕捉するが、特定の構造（例えば骨棘など）または微妙な形状変化を捉えることは不可能である。輪郭全体（全輪郭）を使用すればこの問題は解決するであろう。しかしながら、椎骨の輪郭全体を手作業でアノテートすることは膨大な作業である。従って、自動または半自動で椎骨を領域分割（segment）する多くの方法が提案されてきた。

Gardner等（非特許文献１）は、椎骨の動的輪郭（スネーク（snake））モデリングに基づく半自動システムを開示している。椎骨のデジタル画像上において椎骨の境界線上のポイントがユーザによって指定される。選択された境界線ポイントは、椎骨境界線に自動的にフィットした動的輪郭に対する物理的制約条件として働く。様々なパラメータがモニタされ、スネークに対する適切な自由度を与えるように調整される。自由度が多すぎると、スネークが椎骨の不鮮明なエッジ（辺部）によって誤った方向に導かれる可能性があるので、誤った形状がもたらされることがある。逆に、自由度が少なすぎると、例えば椎骨の角部または骨棘などの、輪郭の鋭い曲がり角が見あたらない状況が起こる可能性がある。

ラテラルＸ線画像における椎骨の半自動領域分割および自動領域分割は、その画像の性質から、困難なタスクである。ラテラルＸ線画像は多くの層（レイヤ）の重ね合わせを示し、このことが関心のある領域（−脊柱（spine）に沿った矢状面−）の区別を難しくしている。これが、多くの古典的な領域分割アプローチ、例えば領域拡張（region growing）に基づくアプローチなど、が失敗してきた理由である。

動的形状モデル（active shape models：ＡＳＭ）は、椎骨の形状および外観に関する事前の知識を利用し、従って画像内の情報にそれほど頼る必要がない、という利点を有する。

Smyth等の論文（非特許文献２）では、椎骨の全輪郭を表すために点分布モデル（point distribution model：ＰＤＭ）が使用され、正常な椎骨と骨折した椎骨を分けるために分類器（classifier）が使用された。６点ランドマーク表現のパフォーマンスと比較して、全輪郭表現のパフォーマンスに改善がみられた。de Bruijne等の論文（非特許文献３）でも同じような結論が下された。

Zamora等（非特許文献４）は、階調レベルエッジ情報を含む動的形状モデル（ＡＳＭ）を使用し、一般化ハフ変換のカスタマイズした実装でそれを初期化した。彼らはその方法をＸ線画像に適用して、腰椎画像の５０％において６．４ｍｍ未満の誤差を実現した。Smyth等も、二重エネルギーＸ線吸収測定法（dual energy x-ray absorptiometry：ＤＸＡ）において椎骨を領域分割するためにＡＳＭ法を使用した。彼らは、ユーザに対して、Ｌ４の最下部、Ｔ１２の最上部およびＴ７の最上部の中間点をアノテートすることを求めている。彼らは健常な椎骨の９５％のケースで１．２３ｍｍ未満のＲＭＳポイントツーライン誤差（root mean square (RMS)point-to-line error）を実現し、骨折の９２％のケースで２．２４ｍｍ未満のＲＭＳポイントツーライン誤差を実現した。

De Bruijne等は形状粒子フィルタリング（shape particle filtering）に基づく完全に自動化された方法を適用し、その結果、平均ポイントツーライン誤差が１．０ｍｍまで低下した。Roberts等（非特許文献５）はＤＸＡ画像内の椎骨を領域分割するために動的外観モデル（active appearance model：ＡＡＭ）と動的オーダリングアルゴリズム（dynamic ordering algorithm）を取り入れた。ユーザにSmyth等と同じ３つのポイントをマークすることを求めることで、彼らは０．７９ｍｍのポイントツーライン誤差を実現した。Roberts等（非特許文献６）によって、より良好な結果が達成されている。彼らは最近、健康な椎骨のレントゲン写真において０．６４ｍｍの平均ポイントツーライン誤差（９５％のポイントが輪郭から２ｍｍ以内にある）を報告した。ユーザは各椎骨のほぼ中央にマークすることが求められている。このアルゴリズムはRoberts等（非特許文献７）によってＤＸＡ画像にまで拡張された。彼らは正常な椎骨において０．６９ｍｍの平均ポイントツーライン誤差を実現し、骨折した椎骨において０．９６ｍｍの平均ポイントツーライン誤差を実現した。

Gardner et al., SPIE Vol. 2710, February 1996, pp 996 - 1008 Smyth et al., Radiology, May 1999, pp 571-578 de Bruijne et al., "Vertebral Fracture Classification," in Medical Imaging: Image Processing, J.P. Pluim and J.M. Reinhardt, eds., Proceedings of SPIE 6512, SPIE Press, 2007 Zamora et al., Medical Imaging: Image Processing, Vol 5032 of Proceedings of SPIE, SPIE Press, 2003, pp 631-642 Roberts et al., Volume 2750 of LNCS, Springer, 2005, pp 733-840 Roberts et al., Proceedings of Medical Image Understanding and Analysis, 2006, Vol.1, pp 120-124 Roberts et al., Proceeding of MICCAI Conference Workshop on joint and bone disease, 2006, Vol. 1, pp 1-8

本発明は、従来技術の上記課題を解決するための手段として、構造物（structure）を含む画像を処理することによって、該画像内の構造物の輪郭（contour）の位置を特定する方法を提供する。本方法は、
前記画像内の構造物はその上に３個乃至１０個のランドマーク（標識点）ポジションがアノテートされており、前記構造物を含む画像を表現するデジタルデータの開始集合（starting set）を用意するステップと、
前記構造物の統計モデルを前記画像上にアノテートされた前記ランドマークポジションにフィットさせ、前記構造物の輪郭の初期推定値を生成するステップと、
推定された輪郭に隣接するポイントから得られた階調レベル情報を使用して、反復的に輪郭境界線を移動させて前記構造物の輪郭の最終推定値を生成するステップと
を含んでなる。

本発明は一般に、画像（例えば、Ｘ線画像、コンピュータ断層撮影（ＣＴ）画像あるいは磁気共鳴（ＭＲ）画像など）内の、一般に予測可能な形状を有する任意の構造物に適用可能である。また画像内においては、画像内の構造物と、画像の背景との間の光強度の違い（コントラスト）が存在する。類似の構造物を含む事前の画像が、後で新たな画像内の構造物の輪郭を位置特定（locate）するために使用される統計モデルのトレーニングに使用できる。

本方法は、身体のパーツや、特に骨（bones）を位置特定するために使用されることがある。この点について、骨は一般に、骨の輪郭の位置を特定するために統計モデルが骨の画像上の手作業によるアノテーションの助けを借りて使用できるような、予測可能な形状を有する。

本方法の好ましい態様において、前記構造物は椎骨（vertebra）であり、前記画像は脊柱（spine）の一部の画像である。本方法は、以下、椎骨の輪郭の位置特定への適用に関連して記述される。しかしながら、本方法は、画像内の任意の構造物の輪郭を抽出すること（−その形状の一般的な予測は、一群の斯かる画像内の類似の構造物の輪郭の統計解析によって推論可能である−）に等しく適用可能であることが理解されるであろう。

例えば、本方法は、手または手首のＸ線画像内の骨の輪郭を位置特定するため、あるいはＭＲもしくはＣＴスキャン画像内の心室もしくは心臓壁の輪郭を位置特定するために使用することも可能である。

椎骨上に複数のランドマークポジション（−統計モデルはこれらのランドマークポジションにフィットさせることができる−）を提供することで、椎骨の縁（エッジ）を位置特定する反復プロセスの好適な開始ポイントが与えられる。好ましくは、椎骨は、その上に４乃至８個、例えば６個のランドマークポジション（−各コーナーに１つ、各椎骨端板の中間点に１つ−）がアノテートされる（打たれる）。椎骨を既に述べたように６個のランドマークでアノテートすることは現行の標準的な実務である。結果として、標準的な実務と比較して、追加の手作業は一切必要とされない。

本方法は、１つの態様として、他の椎骨から成る集合のそれぞれの輪郭を近似するポイントからの情報を用いて椎骨の前記統計モデルをトレーニングするステップを更に含む。椎骨のおおよその輪郭を提供する任意数のポイントが使用できることは理解されよう。しかしながら、好ましい態様では、輪郭を近似するために２０以上のポイント、例えば４０より多くのポイント、５０より多くのポイント、あるいは６０より多くのポイントが使用される。

好ましくは、本方法は、前記他の椎骨から成る集合に属する各椎骨上にアノテートされた、３個乃至１０個のランドマークポジション、例えば４個乃至８個のランドマークポジションからの情報を用いて統計モデルをトレーニングするステップを更に含む。好ましい態様では、他の椎骨から成る集合に属する各椎骨上に６個のランドマークがアノテートされる。

本方法の好ましい態様では、統計モデルをトレーニングする際に使用される椎骨の集合は、未骨折の椎骨と骨折した椎骨を含む。骨折した椎骨をトレーニングに含めることによって、統計モデルは、未骨折の椎骨のほかに骨折した椎骨の縁（エッジ）も位置特定することが可能となる。骨棘を有する椎骨もトレーニング段階に含めることができ、それにより本方法は新たな画像内において骨棘を見つける蓋然性を高めることができる。

本方法の１つの態様として、前記統計モデルは条件付き点分布モデル（conditional point distribution model）である。しかしながら、数学的プロセスが適切に構成されるのであれば、他の統計的形状モデル−例えば、球面調和関数、フーリエ記述子、距離マップ、画像ワーピングおよびvolumetric表現またはmedial表現に基づく統計的形状モデル−が代わりに利用できることは理解されよう。

好ましくは、前記条件付き点分布モデルは、一群の椎骨のそれぞれの輪郭を近似する情報と、前記一群の椎骨に属する各椎骨上にアノテートされた３個乃至１０個、４個乃至８個、または６個のランドマークの情報とから構築される。

追加的に、および／または代替的に、前記条件付き点分布モデルは、前記一群の椎骨のそれぞれの輪郭を近似する情報から構築された第１の点分布モデルと、前記一群の椎骨に属する各椎骨上にアノテートされた３個乃至１０個、４個乃至８個、または６個のランドマークの情報から構築された第２の点分布モデルとから構築される。

本方法の好ましい態様では、前記輪郭の初期推定値は、前記ランドマークポジションにフィットした条件付き点分布モデルの平均値である。

好ましくは、推定された輪郭の反復移動（iterative movement）は、前記輪郭が現在推定値の近傍にあること（近接性）によって制約される。追加的に、および／または代替的に、推定された輪郭の反復移動は、条件付き共分散によって制約される。これらの制約条件によって、初期推定値の周りの探索空間が狭まり、妥当と思われる形状が結果としてもたらされる。

追加的に、および／または代替的に、輪郭境界線の移動は、推定された輪郭に隣接するポイントから得られた階調レベル情報の勾配を、前記統計モデルから得られた等価情報で制約することによって制約される。

好ましくは、輪郭境界線の反復移動は、２回の連続する反復でそれぞれ推定された輪郭の間の差（差分）が事前に設定された限界を下回るまで継続される。例えば、連続する輪郭推定値の間の距離が２ピクセル未満または１ピクセル未満になると、反復プロセスは終了する。

本方法の１つの態様として、各輪郭ポイントを横断する前記輪郭の法線に沿って前記画像内の階調レベル情報をサンプリングすることによって階調レベルプロファイル（grey level profile）が構築される。

本発明は、原則として、デジタル画像から情報を引き出す方法として定義されてきた。しかしながら、当然のことながら本発明は、前記方法を実行するコンピュータのための命令群として、または適切にプログラムされたコンピュータとして、等しく適用可能である。

以下に簡単に説明される添付図面を参照して、本発明の実施形態を以下に詳細に説明する。

６個の初期のランドマークポジションと輪郭がアノテートされた椎骨の一例を示す図である。 本発明の実施の一形態による、椎骨の輪郭の初期推定値の一例を示す図である。 最大許容マハラノビス距離の、結果に及ぼす影響を示す図である。 本発明によって決定された輪郭の実際の輪郭からの距離をヒストグラムと累積分布関数の形で示す図である。 leave-one-out実験からのいくつかの例を示す図である。 ポイント番号に依存する誤差を示す、ポイントインデックスに対する平均ポイントツーライン誤差のグラフを示す図である。 αに対する平方誤差の和のグラフによって、（ａ）は平方誤差の平均和のαに対する依存性を示す図であり、（ｂ）は平方誤差の最大和のαに対する依存性を示す図であり、（ｃ）α＝１．５７における、平方誤差の平均和の、トレーニングセット内の骨折数に対する依存性を示す図である。 標準のＰＣＡとα−ＰＣＡを使用した場合の違いを示す、本発明による領域分割の比較例を示す図である。

以下、本発明を特に脊柱（spine）の椎骨（bertebrae）のＸ線画像の解析に関連して説明する。しかしながら、本発明の方法は、脊柱の他の医用画像−例えば、ＤＸＡ（dual X-ray absorptiometry）、ＣＴ（Computer Tomography）またはＭＲ（Magnetic Resonance）など−に適用できることは理解されよう。

これから述べる全てのステップは医用画像内の任意の構造物の輪郭の抽出に等しく適用可能であり、その形状の一般的な予測は一群の斯かる画像内の類似の構造物の輪郭の統計解析によって推論可能である。例えば、本方法は手または手首のＸ線画像内の骨の輪郭の抽出、あるいはＭＲまたはＣＴスキャンにおける心室や心臓壁の輪郭の抽出に適用されることがある。

本発明の椎骨の輪郭を位置特定するための好ましい方法は２つの主ステップから成る。

一番目のステップは、条件付き点分布モデル（conditional point distribution model：ＰＤＭ）を構築することである。このために、最初に２つの点分布モデル（point distribution model：ＰＤＭ）が構築される。第１のモデルは従来の６つのランドマークポイントでアノテートされた椎骨のトレーニングセットから得られ、第２のモデルは、実際の輪郭アウトラインを近似する、例えば２０個以上の多数のポイントでアノテートされた椎骨の同じトレーニングセットから得られる。次にこれら２つのＰＤＭの関係は、新たなケースに対して、臨床医が手作業でアノテートした６つのポイントのポジションに応じて全輪郭の条件付きＰＤＭを構築することが可能となるようにモデル化される。

次に二番目のステップは、この条件付きＰＤＭを、従来の６つのランドマークポイントでアノテートされた椎骨の新たな画像に適用し、その椎骨の初期輪郭を近似することである。次に、動的形状モデリングを使用し、条件付きＰＤＭ共分散の制約条件に従って初期輪郭を操作して椎骨の実際の輪郭を見つける。

次に、上記２つのステップについて、具体例を用いてより詳しく説明する。

この具体例で使用されるトレーニングセットは、１４２人の患者からの、椎骨の全輪郭と椎骨の６個のランドマークポジションの集合の情報から成る。本例では、画像内の椎骨の撮像と投影の結果として、椎骨は２つの輪郭を持つことが示されたが、下方の輪郭が常に選ばれた。椎骨Ｌ１、Ｌ２、Ｌ３およびＬ４が解析された。従って、計５６８個の椎骨（６４個は骨折した椎骨）がこの調査に含まれた。

画像は１２ビット・ディープで、それらの解像度は５７０ＤＰＩであった。全ての画像はＤＩＣＯＭフォーマットで格納された。アプリケーションはこれほどの高解像度を必要としないので、画像はガウス核で平滑化され、５倍ダウンサンプルされた。

３人の異なる放射線科医師によって６個のランドマークポジションと輪郭がトレーニング目的でマークされた。画像上にランドマークポジションをマークした放射線医師は常に輪郭にもアノテートすることになっていた。６個のランドマークのアノテーションでは、最初にコーナー（角部）がマークされた。続いて、上部コーナーを結ぶ線分の垂直二等分線が表示された。それは放射線医師のためのガイドとしての役割を果たし、放射線医師は最小の高さのポイントにランドマークを配置することになっているが、それがはっきりしない場合は、二等分線のできるだけ近くに配置する。このプロセスは次に下部端板に対して繰り返される。表示された二等分線は、アノテーションプロセス全体にわたって放射線医師同士に食い違いが起こらないように手助けし、ＰＤＭにおける観測者間のばらつきの影響を最小化する。

全輪郭をアノテートするため、放射線科医師は欲しいだけ多数の頂点で多角形ラインを引く。この輪郭は調査のグランドトルース（ground truth）として使用された。図１から分かるように、６個のランドマークと、輪郭が、それより前のアノテーションを示すことなく異なるパス（passes）でアノテートされた。このため、それらは必ずしも重なり合わない。

好ましい実施形態による方法は以下のステップに基づく。
トレーニング：
★放射線専門医によってトレーニング画像上に６個のランドマークと椎骨のおおよその輪郭がアノテートされる。
★データセット内の椎骨がアラインされる。
★２つのＰＤＭが構築される。第１のＰＤＭは一群の椎骨のそれぞれのおおよその輪郭の情報を使用し、第２のＰＤＭは同じ一群の椎骨の各椎骨にアノテートされたランドマークポジションの情報を使用する。

椎骨エッジ（椎骨の縁）の位置特定：
★放射線医師によって椎骨の実際の画像上に６個のランドマークがマークされる。
★画像上の６個のランドマークのポジションと、トレーニング段階で得られた２つのＰＤＭとに基づく条件付きモデルが解析のために構築される。
★条件付きモデルの平均値を用いて、初期解が推定される。
★解を制約するために条件付き共分散（conditional covariance）を使用して、輪郭を画像内の椎骨にフィットさせるために動的形状モデリングが使用される。

次にこれらのステップを更に詳しく説明する。

本発明の実施形態を点分布モデルの利用に関連して説明する。しかしながら、その数学的処理は、他の統計モデル−例えば、球面調和関数、フーリエ記述子、距離マップ、画像ワーピングおよびvolumetric表現またはmedial表現に基づく統計的形状モデル−の利用を可能にするように構成することができることは理解されよう。

点分布モデルを構築する際、輪郭ポイントはトレーニング画像の椎骨上に、椎骨の相当ポイントに配置されなければならない。ポイントの数は放射線医師が任意に選ぶことができるが、しかしながら、ポイントは２０個以上あれば輪郭のアウトラインをマークするのに十分であることが示されている。点分布モデルを構築するとき、輪郭ポイントの数はあらゆるケースで同じである必要がある。従って、放射線医師によってマークされるポイントの数は一貫性を欠く可能性があるので、輪郭においてリサンプリングが必要である。本例で述べたトレーニングセットでは、放射線医師によってアノテートされたポイントの最大数は５３であった。

便宜上、本実施形態では、輪郭モデルは６７個のポイントから成ることが決められた。これにより、ランドマーク間に等しい数のポイントが存在することが可能である。しかしながら、ポイント数が輪郭のアウトラインに十分な精度を確保できるという条件のもとで、異なる、決まった数のポイントが使用できることは理解されよう。従って、６７個のポイントに到達するために、放射線医師がマークしたポイントを使用して輪郭が完成され、６７個のポイントが輪郭に配置された。本実施形態では、初期の６個のランドマークポジションに最も近い輪郭ポイントは、１、１３、２５、４３、５５および６７番目のポイントになるように選ばれた。残りの輪郭ポイントはこれら６個のポイントの間に等間隔に配置された。三番目のセグメント（区分）は、平均的に他の４つのセグメントよりも（およそ）５０％以上長いので、５０％以上多くのポイントを有する。サンプル画像を図１に示す。６個の最初のランドマークポジションは星印でマークされた。また、選択されたポイントによる輪郭はアステリスクでマークされた。

本具体例は６個のランドマークポジションの利用に関して述べられてきた。しかしながら、椎骨の形状を一般になお表現する３個乃至１０個の範囲のランドマークが配置される場合があることは理解されるであろう。

以下のステップにおいて、椎骨の６点表現は、ポイント集合間の並進、スケール（拡大縮小）および回転のずれをなくためにアライン（整列）された。コーナー（角部）からの情報のみを用いて、プロクラステス法（Procrustes method）が使用された。アノテーション形状の複素表現は、次のように定義することができる。

そして、一般化プロクラステス整列（generalised Procrustes alignment）は、次の様に表現することができる。

上式において、各ω_ｉ ^ｐはω_ｉの

上への完全プロクラステス・フィットである。平均形状［μ］の完全プロクラステス推定値は、次の平方積行列（squares and products matrix）の複素和の最大固有値に対応する固有ベクトルとして見いだすことができる。

各椎骨ごとの変換パラメータが全輪郭に適用された。その結果、両方の表現がアラインされた。各形状は次のベクトルで表現された。ただしモデルに応じて、Ｎ＝６またはＮ＝６７である。

トレーニング形状がアラインされたところで、６ランドマークポジション表現と全輪郭表現の両方に対して点分布モデルが構築された。これにより、椎骨Ｌ１、Ｌ２、Ｌ３およびＬ４は単一モデルに混合される。このモデルを構築するため、アライン後のデータベクトルに対して主成分分析ＰＣＡ（Principal Component Analysis）が実行された。ＰＣＡはデータセットの次元数（dimensionality）を減らすために使用することができるテクニックである。それは、データの何らかの射影によって最大分散値が一番目の座標（第１主成分として知られる）に対応し、二番目に大きな分散値が二番目の座標に対応し、・・・（以下同様）の新たな座標系への線形変換である。従って、この新たな表現において第１主成分のみを保持することによってデータセットを単純化することが可能である。この新たな空間の直交基底はデータセットの共分散行列の固有値によって与えられることを示すことができる。

形状モデルは、平均形状

と、保持される固有値λ_ｉ（ｉ＝１，２，．．．，ｋ）および対応する固有ベクトルによって特徴付けられる。また固有ベクトルは、Ｐ^ｔＰ＝Ｉ（Ｐ^ｔはＰの転置行列）を満たすＰ行列（射影行列）の縦列にまとめられる。このとき各形状は次式によって近似することができる。

上式において、ｂは、ある特定の形状に対して計算することができる。

また、誤差ベクトルは次式で与えられる。

ここでｂはｋ個の成分から成る縦ベクトルで、形状の当該モデルの空間への射影を表している。トレーニングセットに沿って、このベクトルの平均値はゼロとなり、その共分散行列Ｃはｋ個の固有値を含む対角行列となる。

ある特定のベクトルｂがもっともらしい形状（plausibel shape）に対応するかどうかを検証するため、それがモデルの平均値、つまりゼロベクトルからあまり離れていないことをチェックしなければならない。それと同時に、もっともらしい形状はゼロベクトルに近いｂベクトルを採用するだけで生成できる。有効領域はｂのマハラノビス（Mahalanobis）距離を制限することによって定義される。限界ｄ_ｍａｘはカイ二乗分布を用いて選ぶことができる。

条件ｄ＜ｄ_ｍａｘが成り立たない場合、ｂベクトルは次の様に修正することができる。

データセット内の骨折（した椎骨の）数が健康な椎骨の数と比べて低いとき、モデルを構築するときにそれらにより高いウェイトを与えることによってモデルへのそれらの影響を増大させた。形状の平均と分散を計算する際に、全体の寄与が等しくなるように、２つの異なるウェイトが正常な椎骨と骨折した椎骨に与えられた。５０４個の健康な椎骨と６４個の骨折した椎骨が利用できたので、ウェイトはそれぞれ

と

であった。

この加重方法の代替として、外れ値（outliers）の重要性を扱う（既に述べた）ＰＣＡの変形態様を導入することができる。小さな部分集合のサンプルが外れ値として現れることがあり、信頼できそうなデータになお相当すると考えられるケースが存在する。これは、椎骨の得られたデータの一部が、骨折した椎骨を含む場合に当てはまる。これらは外れ値として現れることがあるが、それらが提供する情報はなお重要である。このケースでは、外れ値のモデリングは重要であり、むしろそれらを無視せず、結果を向上させるためにそれらを使用することができる。

ＰＣＡは最小自乗の意味でデータを最適に近似する部分空間を張る直交線形変換である。これは変換後の座標の分散を最大化することによって遂行される。

データの次元数をＮまで減らさなければならない場合、ＰＣＡの等価形式は、データセットにおける元のデータポイントを再現する際に生じる誤差を最小化するＮ−セットの直交ベクトル（これらはＰ行列にまとめられる）を見つけることである。その誤差は次の様にＬ_２ノルムで測られる。

上式において、Ｃはコスト関数、Ｎはトレーニングケースの数、ｘ_ｉは近似するための中心化されたデータベクトル（centred data vectors）である。

最小自乗法は、データセットの中に外れ値が存在する場合、それらが望ましいソリューションからの結果を歪め、統計推定値の誤り率と歪みを膨らませる可能性があるので、ローバスト性が低い（非特許文献［Hampel et al. : “Robust statistics: the approach based on influence functions" Wiley 1986］参照）。このため例えば、上式のコスト関数を、例えばデータサンプルが外れ値と見なされるときにゼロとなるバイナリ変数を導入することによって修正することによって、ＰＣＡに対する外れ値の影響を低減する試みが行われてきた。

上式において、Ｖ_ｉはバイナリ変数の集合である。項η（１−Ｖ_ｉ）は、最適化によって全てのｉに対してＶ_ｉ＝０となる自明解に収束することを防止する。

これとは対照的に、かつ、通常のＰＣＡにおけるように平方データ再現誤差を直に最小化することとは対照的に、提示するΦ−ＰＣＡアルゴリズムは次式を最小化する。

上式において、ΦはΦ（ｘ^２）が凸型になるような二回微分可能関数である。Φが二回微分可能であるという事実から、最適化において、二次収束をもたらすヘシアン（Hessian）ベースの方法を用いることが可能である。凸性要件はＣに対し、ただ一つの最小値が存在することを保証する。

シンプルかつ同時に強力な関数形は、

である。ただし、凸性条件を満足するにはα＞０．５でなければならない。この特殊なケースをα−ＰＣＡと呼ぶことにする。αのより大きな値（一般にα＞１）は、通常のケースと比べてよりコスト高であるので、外れ値を引き立たせることになる。特に、α＝∞は

ノルムの最小化をもたらし、このため、Ｌ_２ノルムで測られた、形状に関する最大の再現誤差をもたらす。他方、より小さな値（０．５＜α＜１）は、よりローバスト性の高いＰＣＡをもたらすという正反対の効果を生む。α＝０．５の場合はＬ_１ノルムを最小化する。最後に、α＝１は標準のＰＣＡに相当する。

データポイントｘ_ｉは中心化されなければならない。このことはそれらの平均値がそれらから差し引かれなければならないこと、つまりｘ_ｉ＝ｓ_ｉ−μを意味する。ここでｓ_ｉは元々の非ゼロ平均のデータサンプルを表す。本提案のアルゴリズムにおいて、“平均”は標準のＰＣＡにおけるようなデータポイントの成分単位の算術平均ではもはやなく、次式を最小化するベクトルである（データポイントをＭ次元と仮定する）。

ベクトルｘ_ｉの計算が完了したら、当初の式におけるコスト関数を最小化する基底ベクトル（を縦列に持つ直交行列）Ｐを探索することから成るφ−ＰＣＡが実行できる。

またはＰの閉形式表現が存在しないために、ＣとＣ_μの両方を最小化するには、数値的方法が必要とされる。

標準的なＰＣＡとφ−ＰＣＡとの重要な違いは、後者において、所望の次元数が変更される場合に主成分が再計算されなければならない点である。他方、標準的なＰＣＡでは、Ｎ_２＞Ｎ_１を仮定すると、斯かるＮ_２とＮ_１に対する解析において最初のＮ_１個の主成分は共通している。

コスト関数

の勾配（gradient）とヘシアン（Hessian）の表式は全くシンプルで、すぐに計算できる。成分単位の算術平均を初期設定として使用すると、ニュートン法は次の解にすぐに収束する。

上式において、勾配

は次式のように一次偏微分から成る列ベクトルである。

また、ヘシアン行列Ｈは、次式のように二次偏微分から成る。

平均値がデータポイントから差し引かれた後、当初のコスト関数Ｃが最小化されなければならない。この関数は、Ｐ^ｔＰ＝１を満足するような直交行列Ｐに対してグローバルミニマム（global minimum）に近づくという、興味深い特性を有する。これにより、最適化の際にＰをこの条件を満足するように制約する必要がなくなることが可能となる。このことは先ほど言及した特性のおかげであって、一般には、この条件を置かないということは、Ｐ^ｔＰが射影行列を表さず、このためＰＰ^ｔｘ_ｉ−ｘはもはや再現誤差を表さないことを暗に意味する。

最小化問題では、勾配の表式のみが導入される。というのは、ヘシアン行列の表式はあまりにも複雑で、その計算は高くつき過ぎるからである。行列計算を使用すれば、全ての偏微分をいっせいに計算することが可能である。

勾配がわかると、異なる標準的技術を使用してＰを更新することができる。シンプルな勾配降下スキーム（gradient descent scheme）では、例えば、次の様な式が用いられる。

上式において、ｋはステップサイズ（刻み幅）である。

通常のＰＣＡでは、最適なＰを迅速に確保するために、ラインサーチ（line search）が初期設定として使用できる。このアルゴリズムでは、各反復ごとに異なるステップサイズが調べられ、コスト関数Ｃの最小値をもたらすステップサイズが確保される。直交性条件は、コスト関数と勾配の表式を単純化するものではあるが、Ｐ行列は（直交行列に収束するとはいえ）無制約に変更されているので、プロセス全体を通して仮定することができないことに言及することは重要である。

形状モデル（Cootes et al. 1995）では、一群のランドマークが、一群の既にアライン済みの形状の上に定義される。ランドマークのｘ座標とｙ座標を束にすることによって、１つのデータベクトルｓ_ｉが形状ごとに構築される。次に、ある程度の近似誤差を代償として、イクスプリシットなデカルト座標よりも、より次元数が低く、より特異性が高い表現で形状を表現することを目的として、平均値がデータベクトルｓ_ｉから差し引かれ、結果のｘ_ｉデータベクトルにＰＣＡが実行される。標準のＰＣＡとφ−ＰＣＡに基づく形状モデルの違いについて説明する。

最初に、プロクラステス法を用いて形状がアライン（整列）され、それらの平均値が計算される。形状をアラインするために、回転（rotation）、並行移動（translation）および拡大縮小（scaling）が許される。次式を最小化することによって、整列パラメータと平均値が同時に最適化される。

上式においてＴ_ｉ（ｚ_ｉ，θ_ｉ）は一群のパラメータθ_ｉに基づくアライン後のｓ_ｉ形状を表す。制約条件μ^ｔμ＝１は形状がゼロに収縮するのを防止する。Cootes等の文献に記述された反復アルゴリズムが以下の問題を解決するために使用された。
１．第１の形状のサイズを規格化し、それを平均の第１推定値として使用する。
２．全ての形状を平均の現在推定値にアラインする。
３．アライン後の形状の平均値を見つけることによって、平均の推定値を更新する。
４．平均の新たな推定値のサイズを規格化する。
５．収束するまでステップ２に進む。

第３ステップにおける平均値はコスト関数を数値的に最小化することによって見つけなければならない。しかしながら、φ（ｔ^２）の最小化はｔ^２＝｜｜ＰＰ^ｔｘ−ｘ｜｜^２の最小化と同等であるので、第２ステップにおけるアラインメント（整列）は、標準的な方法で自乗距離の和を最小化することによって、容易に計算することができる。この特徴のもう１つの結果は、ある形状のｂ_ｉのＰＣＡ座標は、次の様に、通常のＰＣＡにおけるものと同じ方法でなお計算することができるということである。

全輪郭モデルＦの主成分の分布は、ランドマークポジションＬの主成分に応じて、条件付きガウス分布としてモデル化することが可能である。μ_Ｆ、μ_Ｌ、Σ_Ｆ、Σ_Ｌが、トレーニングデータ全体にわたる、２つのモデルの主成分の平均値と共分散行列であり、Σ_ＦＬとΣ_ＬＦがそれらの相互共分散行列である場合には、次の様な式を書き下すことが可能である。

上式において、６個のポイントの主成分座標Ｌが与えられた全輪郭の主成分座標に対する、μは条件付き平均、Σは条件付き共分散行列である。

主成分の代わりにポイントとランドマークのポジションをモデル化することも可能である。後者は共分散行列のサイズがより小さいときに利点がある。ポイントの座標が直に使用される場合、Σ_ＬＬはポジションに関する多重共線形性のせいで非可逆になる可能性があり、正則化が必要になるだろう。

図２に示す様に、条件付きモデルの平均値は初期設定として使用することができ、共分散はモデルを画像にフィットさせるときに有用である。条件付き共分散は一般にＣの無条件共分散よりずっと“小さい”。分布の微分エントロピーは無条件モデルから条件付きモデルまで対数単位でほぼ１０ポイント減少する（−１３．８８から−２３．８５まで）。従って、新たな条件付き共分散に基づいて定義されたある特定の値のマハラノビス距離によって制限された領域内において、条件付き平均の周りで解を探すことが可能である。そのために探索空間は狭まり、そのおかげでモデルをフィットさせることがより容易になり、かつ、形状が６個のランドマークポジションから比較的遠く離れる可能性が低くなる。

６個より多くのランドマークまたは６個より少ないランドマークを有する実施形態では、これらのランドマークポジションは、条件付き共分散を計算するときに、本実施形態におけるものと同じ方法で使用される。

本例では、トレーニングデータ内の６個のランドマークは輪郭上のそれらの対応するポイントに留まることは制約されない。この制約がないことから、１１−Ｄ全輪郭条件付き形状モデルは、マハラノビス距離がより長く、従って蓋然性がより小さくはなるが、条件付きでないモデルと同じ形状を正確に表すことが可能である。このことは重要である。なぜならば、それにより条件付き共分散行列がフルランク、つまり可逆的であり続けることが可能となるからである。

動的形状モデル（active shape model：ＡＳＭ）は、形状モデルを画像内の椎骨の輪郭にフィットさせようと試みる反復アルゴリズムである。最初のステップは、６個の与えられたランドマークポジションのコーナーを形状モデルの平均にベストフィットさせる、並進（ｔ_ｘ，ｔ_ｙ）、回転（θ）およびスケール（ｓ）パラメータを見つけることである。これらの姿勢（pose）パラメータは、画像内のポイントの（“物理的”座標における）ポジションＸと、形状モデルの“規格化”座標におけるそれらのポジションｘとの間の切り替えを可能にする変換を定義する。姿勢パラメータはプロセス全体を通して一定に保たれる。

“規格化”座標において６個のランドマークポジションで計算される初期解から出発して、輪郭の法線に沿った並進が各反復ごとにモデル内のあらゆる輪郭ポイントに対して提案される。それを計算するため、一群の候補ポジションｔが各輪郭ポイントにおいて輪郭の法線に沿って選択される。各ｔごとに濃淡値プロファイルが、画像内の階調レベルを輪郭の法線に沿ってｔの周りでサンプリングすることによって構築される。プロファイル内のポイントに対して微分（derivative）が計算され、次いで、微分プロファイルの絶対値の和が１になるように拡大縮小（スケール）される。これによりアルゴリズムはコントラスト変動に対して頑強（ローバスト）なものとなる。次に結果のプロファイルｐ_ｉ（ｔ）が、輪郭上の同じ輪郭ポジションにおける、トレーニングケースのものと比較される。

ｐ_ｉ（ｔ）が、輪郭ポイントｉの周りの、プロファイルにおけるポイントｔの周りの（セミレングスＴ_ｔを用いて、例えば区間［ｔ−Ｔ_ｔ，ｔ＋Ｔ_ｔ］における）規格化された微分（derivatives）のベクトルを表すとすれば、ｐ_ｉ（ｔ）をトレーニング例から構築されたモデルと比較することによって、各ｔごとにフィットネス関数ｆ_ｉ（ｔ）が計算できる。

上式において、

はトレーニングケースにおけるポイントｉの周りの長さ２Ｔ_ｔ＋１のプロファイルの平均値である。またＳ_ｐｉはプロファイルにおける各要素の（独立な）分散を含む対角行列である。関数ｆ_ｉ（ｔ）はちょうど区間［−Ｔ_ｐ，Ｔ_ｐ］で定義される。このため、Ｔ_ｐは、輪郭の探索が現在のポイント位置からどれだけ離れて実行されるかを定義する。ｔ_ｍがｆ（ｔ）を最小化する場合、アルゴリズムがよりスムーズな仕方で発展するようにシフト（２／３）ｔ_ｍが選ばれる。

この新しい所望の形状Ｘ＋ｄＸは規格化座標に変換され、ｘ＋ｄｘになる。形状モデルパラメータｂは可能な限りｘ＋ｄｘにフィットするように更新される。

上式において、Ｗ_ｓはフィッティングにおける各ポイントの重要度を測るウェイトの対角行列である。ウェイトは変位の大きさとフィットの良さに依存する。

輪郭を更新する前に、ｄｂが妥当と思われる形状をもたらすことをチェックすることが重要である。条件付きモデルの共分散における情報が有用となるのはこの時点である。

そして、

従って、ｄ_ｍａｘはアルゴリズムがどれだけ自由に輪郭を画像内のエッジにフィットさせられるかを制御するパラメータである。値が大きいと、結果は主成分空間を動き回ることができ、画像内おいてエッジが不鮮明な場合には、あり得ない解がもたらされる可能性がある。値が小さいと、アルゴリズムはほとんどモデルに依存することになり、分布の平均により近い、より保守的な結果がもたらされる。これは、アルゴリズムが、特に骨折または骨棘（こっきょく）を含む異常なケース−図３に示す様に、解は主成分空間において初期設定から比較的離れている−において、正しい解を見つけることを妨げる可能性がある。具体的には、図３は、最大許容マハラノビス距離の、結果に及ぼす影響を示している。図３（ａ）では、形状モデルは輪郭を骨棘にフィットさせることは不可能である。図３（ｂ）では、閾値は１．５だけ増大されており、輪郭は骨棘をより良く近似する。

新座標は次式を用いて容易に計算される。

そして、この新座標は、−ｔ_ｘ、−ｔ_ｙ、−θ、ｓ^−１によって定義される変換によって、物理的座標に逆変換される。次いで、新たな反復を開始して、各ポイントごとにもう一度、新しい並進が提案される。２つの連続する反復における形状の差が、ある特定の限界を下回る場合には、プロセスは終了する。例えば、連続する推定された輪郭の間の距離が２ピクセルないしは１ピクセル未満の場合、反復プロセスは終了する。

まず最初に、アルゴリズムにおける異なるパラメータに対する適切な値を見つけるために、いくつかの予備的な実験が行われた。保持すべき主成分の数の選択に関して言えば、それは維持すべき分散の割合に依存する。６個のランドマークと全輪郭に対して、それぞれ７個の成分と１１個の成分があれば、全分散のほぼ９７％をキープするのに十分であった。そこから割合を引き上げることは、主成分の数の点で非常に高くつくものとなる。

プロファイル長に関しては、濃淡値モデルのプロファイルセミレングスＴ_ｔ＝１２ピクセル（２．７ｍｍ）を使用し、Ｔ_ｐ＝８ピクセル（１．８ｍｍ）とすることは良い結果をもたらすことが見いだされた。Ｔ_ｐは現在の解からどれほど離れたところで輪郭を見つけようとするかを表す。このパラメータを大きくし過ぎると、探索プロファイルは長くなり過ぎることになり、その結果、アルゴリズムが間違ったエッジを捕捉する蓋然性が、特にこのエッジがあり得ない形状を表さない場合に、より大きくなる。これは一般的に、図５（ｂ）に示したタイプの“二重輪郭（double contour）”ケースで起こる。

最後に、最大マハラノビス距離ｄ_ｍａｘの選択が自由度１１に対するχ^２分布の助けを借りて行われた。１０％のテール確率を任意に固定して、４．１に等しい妥協リミットが選択された。特に断りがない限り、全てのパラメータは全ての更なる実験で固定した状態に保たれた。

leave-one-out実験が実施された。各椎骨ごとに、その他の画像からコンプリートモデル（６点ランドマーク、全輪郭、プロファイル）が構築され、輪郭の位置特定に使用された。医師によってアノテートされた輪郭の各ポイントと検出された輪郭における最も近いポイントとの間の距離（ポイントツーライン距離（point-to-line distance））が各椎骨ごとに計算された。図４に誤差の分布を示す。具体的に、図４は実際の輪郭までの距離をヒストグラムと累積分布関数の形で示している。

ＲＭＳ（二乗平均平方根）誤差は０．６８ｍｍに等しく、それに対して平均誤差は０．４８ｍｍであった。またその標準偏差は０．４８ｍｍであった。手作業でアノテートされた輪郭の１ｍｍ以内にポイントの８９％が存在、１．５ｍｍ以内にはポイントの９６％が存在、２ｍｍ以内にはポイントの９８％が存在していた。各椎骨における最大誤差の平均は１．５３ｍｍであった。

モデルを構築する際に、骨折に、より大きなウェイトが与えられたという事実によって、骨折した椎骨に対応する形状にフィットさせるときにモデルのパフォーマンスが低下し過ぎることが防止される。骨折に対する平均誤差は０．５４ｍｍであった。ポイントの８６％は骨折における実際の輪郭の１ｍｍ以内に存在し、１．５ｍｍ以内に９４％が存在した。

結果は、特に骨棘および骨折を含む、輪郭が平滑化し過ぎる可能性のあるケースで、視覚的にも診断された。図５にいくつかのサンプル画像を示す。図から分かるように、モデルは、グランドトルース（ground truth）が（骨棘の周りのように）急カーブを描くときに、そのようなグランドトルースを近似する難しさを持つが、それ以外の場合には、定量的な結果が示唆するように、全くうまく機能する。

結果は異なる調整後パラメータの変化にあまり敏感でないことに注意することも重要である。主なパラメータを変化させるときの平均誤差のインクリメントは、
＊モデルの構築に使用されるプロファイル長に対しては：
選ばれた値の周りの４ピクセル幅領域において４μｍである。
＊輪郭探索で使用されるプロファイル長に対しては：
選ばれた値の周りの４ピクセル幅領域において６μｍである。
＊カットオフされたマハラノビス距離に対しては：
選ばれた値の周りの１ユニット幅領域において２５μｍである。

これらの値は本方法がパラメータの非最適な選択に対してローバスト性を持つことを示唆している。カットオフされたマハラノビス距離は、結果に最も影響を及ぼすパラメータであり、自由度（外れ値に対するより良い近似）と安全性（妥当と思えない形状のより低い蓋然性）の間のトレードオフを表す。プロファイル長は主に収束速度に影響を与える。

好ましい実施形態では、モデルを構築するときに骨折に対してより高いウェイトが与えられる。全体の平均パフォーマンスは低下するが、最も難しいケースにおける精度は向上し、よって最大誤差が低下する。それにもかかわらず、０．５ｍｍ未満の平均ポイントツーライン誤差が実現された。誤差の９５％は１．３６ｍｍ未満であった。パフォーマンスは既に述べた最も最近のシステムのパフォーマンスよりも高い。

６個のポイントと全輪郭（full contour）を異なるパス（passes）でアノテートすることによって、解は６個のランドマークを通ることを制約されない。これにより形状モデルにフレキシビリティが加わり、その結果、その形状モデルがより広範な異なる形状にフィットすることが可能となる。異なる医師がアノテーションを行った事実もモデルにある興味深い可変性を加える。

図６はポイント番号に依存する誤差を示している。６個のランドマークに対応するポイントは星印★でマークされ、手作業で置かれたランドマークから本当の輪郭までの距離は×印でマークされている。ここで、平均ポジション誤差は、骨棘の一般的な位置である、三番目のランドマークの直後、かつ四番目のランドマークの直前に明らかなピークを有することは明らかである。曲線は、下部端板の中間点を除いて、ランドマークに対応するポイントの周りに極小点を有する。６個のランドマークは輪郭上に拘束されていないので、これは可能である。これらのポイントから輪郭までの平均距離は同じ図面において×印でマークされている。

実施の一形態において、骨棘を含む更なるトレーニングケースが使用できる、あるいは、骨折と同様に、より高いウェイトを斯かるケースに与えることができる。代わりに、より高いフレキシビリティがＡＳＭに対してポイント２７とポイント４１（骨棘の一般的な位置）の周りで、例えば輪郭をある種のスムーズな曲線でポイントにフィットするように修正することによって許されることがある。より複雑なソリューションは骨棘の周りのポイントの対応を改善する可能性がある。所与のデータセットに対するベストな選択基準（hypothesis）は最大の圧縮をもたらすものであるという最小記述長（minimum description length：ＭＤＬ）が、対応するポイントを選ぶことに使用できる。

代わりのデータセットの結果もα−ＰＣＡを使用して外れ値に適用される代わりのウェイティングに基づいて与えられる。調査は１４１人の患者の脊柱からのラテラルＸ線から成るデータセットに基づく。椎骨Ｌ１乃至Ｌ４が３人の異なる放射線専門医によって輪郭が描かれ、調査のground truthが提供された。ＭＤＬアルゴリズムを用いて各椎骨ごとに６５個のランドマークが抽出された（非特許文献［Thodberg “Description Length Shape and Appearance Models", proceedings of Information Processing in Medical Imaging (2003) Springer］参照）。同じ放射線医師は椎骨の骨折タイプ（wedge, biconcave, crush）と等級（軽度、中度、重度）に関する情報も提供した。加えて、彼らは標準的な６点形態計測で使用される６個のランドマーク（−コーナー（角部）と、両方の椎骨端板の中間点とに位置する−）もアノテートした。これらのポイントは前高（anterior height）、中高（middle height）、後高（posterior height）を定義し、これらは骨折の等級とタイプを推定するために使用される。

通常のＰＣＡと（異なる値のαに対する）α−ＰＣＡは両方とも、全てのケースでデータの全分散のほぼ９５％を維持することが可能な７つの（α−）ＰＣＡ座標をキープするデータセットに適用された。両方のアルゴリズムに対して、平均と最大の平方再現誤差も計算された。誤差の、トレーニングセットにおける骨折数に対する依存度も調べられた。成分数が増せば、より良い精度が実現されるとともに、モデルの特異性に関する良好なトレードオフがなお実現されるが、本実験では、ＰＣＡとα−ＰＣＡとの間の違いをより良く示すために、より小さな成分数を保ったことに留意すべきである。

最後に、画像内の椎骨Ｌ１−Ｌ４を領域分割するための動的形状モデルでＰＣＡとα−ＰＣＡがテストされた。２つの形状モデルが構築された。１つは６個のランドマークに対するもので、もう１つは全輪郭に対するものである。また両モデルの（α−）ＰＣＡ座標の間の関係は条件付きガウス分布にフィットさせた。モデルにおけるより高いフレキシビリティを可能にするため、より大きな主成分数が利用された。６個のランドマークに対する主成分数は７、全輪郭に対する主成分数は１１とした。これらは両方のケースで全分散の約９８％をキープする。

条件付き分布の平均値が初期設定として全輪郭の領域分割（segmentation）に使用された。各反復ごとに、以下のルーピングにおいて各ポイントごとの所望のポジションを計算するために、輪郭に垂直なプロファイルに沿った階調レベル情報が使用された。モデルを新たなポイントにフィットさせることによって新たな輪郭が計算できる。新たな（α−）ＰＣＡ座標ｂから条件付き平均までのマハラノビス距離を測るために条件付き共分散が使用された。それがある特定の閾値Ｄ_ｍａｘを超える場合には、Ｄ（ｂ）≦Ｄ_ｍａｘを保証するために、ベクトルがスケールダウンされたｂ’＝ｂ（Ｄ_ｍａｘ／Ｄ（ｂ））。このようにして、解は６個のランドマーク近くに留まるよう制約された。このプロセスは収束するまで繰り返される。

図７（ａ）と図７（ｂ）は、ラベルされたポイントにモデルをフィットさせるときの、平方誤差の和（sum of squared errors：ＳＳＥ）のαに対する依存性を示している。予測されるように最大誤差はαの増加とともに減少し、骨折に対してはそれは更に加速する。平均誤差は１より小さなαの値が、骨折の場合における誤差を、それらはもはやモデルにおいて重要ではないので、どのように増大させていくかを示すとともに、未骨折の場合においては、それほどではないにしても、それをどのように減少させていくかを示している。未骨折の椎骨は一般に既に全くうまくモデル化されている。１より大きな値は骨折における結果を、未骨折の椎骨においてそれらを若干悪くはするが、最初に改善する。最後に、αが増大しすぎる場合、モデルは最も起こりそうにないケースにしかフィットしない傾向があり、平均の結果は未骨折の椎骨と軽度の骨折の両方に対して悪化する。

本実験では、未骨折の椎骨と、利用可能な骨折の総量の異なる割合（１２．５％乃至１００％で１２．５％のインクリメント）の全てでモデルが構築された。αは１．７５に設定された。これにより、既に示した結果によれば、骨折した椎骨と未骨折の椎骨における最大誤差と平均誤差の良好なトレードオフが実現された。図７（ｃ）は、トレーニングセット内の骨折数が比較的小さいときには、α−ＰＣＡが通常のＰＣＡを明らかに凌ぎ、特に有用であることを示している。

椎骨Ｌ１乃至Ｌ４は、放射線医師がアノテートした６個のランドマークで調整された形状モデルを使用する利用可能な画像から、標準のＰＣＡとα−ＰＣＡ（α＝１．７５）の両方を用いて領域分割された。実験は、ある特定の画像の領域分割に使用されるモデルが全てのその他のモデルに基づいて構築されるleave-one-out方式で実施された。

本当の輪郭からアルゴリズムの出力までのポイントツーライン誤差を以下の表１に、p t-test（正確にはpaired, double-sided t-test）とp signed rank（正確にはpaired, double-sided Wilcoxon signed rank test）からの結果であるｐ値（ｐ-values）と一緒に示す。結果は、標準のＰＣＡは未骨折の椎骨のケースにおいてより低い平均誤差をもたらすが、α−ＰＣＡは、全平均誤差が若干増大するものの、骨折の重症度の異なる等級に沿ってより均一な結果を与えることを示している。さらに、α−ＰＣＡは骨折のケースにおいて標準のＰＣＡを、特に重症な場合に、大きく凌ぐ。それは、トレーニングデータに骨折情報を必要とすることなく、異なる重要度を各ケースに連続式に割り当てる特徴も有する。仮にこの情報が利用可能であったとしたならば、２つの異なるモデルを構築することが可能であったであろう。しかし、その場合、多数のトレーニング用の骨折が必要とされたであろう。ｐ値（ｐ-values）に関しては、両方のテストは、２つのセットアップの間の平均値の違いが有意であることを示している。

最後に、図８には、標準のＰＣＡとα−ＰＣＡ（α＝１．７５）で領域分割された２枚のレントゲン写真が、放射線医師によって与えられた輪郭と一緒に表示されている。それらは共に重症な骨折に相当する。特に形状がその方向を急激に変えるポイントの周りで、α−ＰＣＡは実際の形状のより良好な近似を実現する。

上述した実施形態の変更および変形は、本願発明の範囲内で行えることが理解されるであろう。

Claims (16)

1. 構造物を含む画像を処理することによって、該画像内の構造物の輪郭の位置を特定する方法であって、
   前記画像内の構造物はその上に３個乃至１０個のランドマークポジションがアノテートされており、前記構造物を含む画像を表現するデジタルデータの開始集合を用意するステップと、
   前記構造物の統計モデルを前記画像上にアノテートされた前記ランドマークポジションにフィットさせ、前記構造物の輪郭の初期推定値を生成するステップと、
   推定された輪郭に隣接するポイントから得られた階調レベル情報を使用して、反復的に輪郭境界線を移動させて前記構造物の輪郭の最終推定値を生成するステップと、
   を含んでなる、半自動式輪郭検出方法。
2. 前記構造物は骨であることを特徴とする請求項１に記載の半自動式輪郭検出方法。
3. 前記構造物は椎骨であり、前記画像は前記椎骨を含む脊柱の一部分の画像であることを特徴とする請求項１または２に記載の半自動式輪郭検出方法。
4. 他の椎骨から成る集合のそれぞれの輪郭を近似するポイントからの情報を用いて椎骨の前記統計モデルをトレーニングするステップを更に含むことを特徴とする請求項３に記載の半自動式輪郭検出方法。
5. 前記他の椎骨から成る集合に属する各椎骨上にアノテートされた３個乃至１０個のランドマークポジションからの情報を用いて前記統計モデルをトレーニングするステップを更に含むことを特徴とする請求項４に記載の半自動式輪郭検出方法。
6. 前記統計モデルのトレーニングで使用される前記椎骨集合は未骨折の椎骨と骨折した椎骨を含むことを特徴とする請求項４または５に記載の半自動式輪郭検出方法。
7. 前記統計モデルは条件付き点分布モデルであることを特徴とする請求項１乃至６のいずれか１項に記載の半自動式輪郭検出方法。
8. 前記条件付き点分布モデルは、一群の椎骨のそれぞれの輪郭を近似する情報と、前記一群の椎骨の各椎骨上にアノテートされた３個乃至１０個のランドマークの情報とから構築されることを特徴とする請求項７に記載の半自動式輪郭検出方法。
9. 前記条件付き点分布モデルは、一群の椎骨のそれぞれの輪郭を近似する情報から構築された第１の点分布モデルと、前記一群の椎骨の各椎骨上にアノテートされた３個乃至１０個のランドマークの情報から構築された第２の点分布モデルとから構築されることを特徴とする請求項７または８に記載の半自動式輪郭検出方法。
10. 前記条件付き点分布モデルは、処理されている画像内の前記ランドマークポジションの位置に依存する条件付きガウス分布であることを特徴とする請求項７乃至９のいずれか１項に記載の半自動式輪郭検出方法。
11. 前記条件付き点分布モデルは、処理されている画像内の６個のランドマークポジションの主成分座標に応じて、一群の椎骨のそれぞれの輪郭を近似する情報から構築される点分布モデルの主成分をモデリングする条件付きガウス分布であることを特徴とする請求項１０に記載の半自動式輪郭検出方法。
12. 前記輪郭の初期推定値は、前記ランドマークポジションにフィットした条件付き点分布モデルの平均値であることを特徴とする請求項７乃至１１のいずれか１項に記載の半自動式輪郭検出方法。
13. 推定された輪郭の反復移動は、条件付き共分散と、条件付き平均が前記輪郭の現在推定値の近傍にあることによって制約されることを特徴とする請求項１乃至１２のいずれか１項に記載の半自動式輪郭検出方法。
14. 前記輪郭境界線の移動は、推定された輪郭に隣接するポイントから得られた階調レベル情報の勾配を、前記統計モデルから得られた等価情報で制約することによって制約されることを特徴とする請求項１乃至１２のいずれか１項に記載の半自動式輪郭検出方法。
15. 前記輪郭境界線の反復移動は、２回の連続する反復でそれぞれ推定された輪郭の間の差が事前に設定された限界を下回るまで継続されることを特徴とする請求項１乃至１４のいずれか１項に記載の半自動式輪郭検出方法。
16. 階調レベルプロファイルが、各輪郭ポイントを横断する前記輪郭の法線に沿って前記画像内の階調レベル情報をサンプリングすることによって構築されることを特徴とする請求項１乃至１５のいずれか１項に記載の半自動式輪郭検出方法。

Images