JP2007203040A - コンピュータ支援小結節検出のための局所肺構造物分類のシステム及び方法 - Google Patents

Abstract

【課題】局所肺構造物分類のシステム及び方法において、高い融通性と低い計算コストでもって偽陽性の数を減らすこと。
【解決手段】デジタル画像において肺構造物を分類する方法において、３次元（３Ｄ）デジタル画像内の１つ又は複数の目標構造物の近似的な目標構造物位置を提供するステップと、前記の近似的な目標構造物位置に関して異方性ガウシアンモデルを当てはめ、より詳細な３Ｄ目標モデルと１つ又は複数の目標構造物の中心位置とを生成するステップと、３Ｄ目標モデルの各々を３Ｄ球体へと変形させるステップと、変形された３Ｄ球体の各々に関して境界多様体を構成するステップと、境界多様体上でクラスタを識別し、前記１つ又は複数の目標構造物を分類するステップを行う。
【選択図】図７

Description

本出願は２００６年１月２５日提出のＢａｈｌｍａｎｎ他による米国仮出願第６０／７６１，９２７号"Local Pulmonary Structure Classification for Computer-Aided Nodule Detection"の優先権を主張するものであり、前記米国出願の内容は参照により本願の開示内容に含まれるものとする。
本発明はデジタル医用画像における局所構造物の種類の分類に関するものである。

昨年、アメリカ合衆国において、１６０，０００人を超す人々が肺癌で亡くなっている。肺癌を防ぐには禁煙が最良であるが、患者の診断を改善するには早期発見が鍵である。癌が早期に発見され、手術が行われれば、Ｉ期の非小細胞肺癌の患者の５年生存率は６０％〜８０％である。しかしながら、手術を受けない患者は５年生存率が１０％でしかない。

コンピュータ断層撮影（ＣＴ）スキャンのようなイメージング技術は早期発見のための非侵襲的かつ高感度の手法を提供する。胸部ＣＴスキャンにおける肺小結節コンピュータ支援検出及び診断（ＣＡＤ）は人為的ミスの可能性を低減し、より効果的で標準化された診断プロセスを提供する。ＣＴスキャンでは、肺小結節は様々な形と大きさの稠密な塊として現れる。これらは血管や胸膜のような他の構造物から孤立している場合もあれば、付着している場合もある。

最近、薄スライスＣＴにおける小結節の自動検出及び分類に関して、領域成長及び自動閾決定や、小結節を選択し雑音を抑制する３Ｄフィルタ、小結節のマッチング、及び、変形可能な幾何学的輝度テンプレートを使用したガウシアン小結節モデルによるテンプレートマッチングを含む多くの技法が提案されている。しかしながら、これら従来技術によるＣＡＤシステムの欠点は小結節と血管のような他の稠密な構造物との区別にある。ほとんどのシステムにおいて円形性の仮定が採用されているため、湾曲した血管やそれらの接合部がしばしば誤って小結節として検出され、結果として偽陽性（ＦＰ）の症例となる。

このような偽陽性の数を減らすため、すでに２種類の解決法が提案されている。１つは、血管のツリー再構成のための曖昧形状分析を用いた相関をベースとしたフィルタにより関心領域を強調することであり、もう１つは、ヘッシアンに基づいた分析により得られる血管中心軸のトラッキングを利用することである。前者のアプローチには融通性の欠如という欠点がある。研究室で用いられる単純な構造物テンプレートでは多くの複雑な血管形状や血管トポロジーを処理することができないだろう。他方で、後者のアプローチはより不規則な構造物を処理することができるものの、計算量の点で非常にコストがかかる。

本発明の課題は、局所肺構造物分類のシステム及び方法において、高い融通性と低い計算コストでもって偽陽性の数を減らすことである。

上記課題は、デジタル画像において肺構造物を分類する方法において、３次元（３Ｄ）デジタル画像内の１つ又は複数の目標構造物の近似的な目標構造物位置を提供するステップと、前記の近似的な目標構造物位置に関して異方性ガウシアンモデルを当てはめ、より詳細な３Ｄ目標モデルと前記１つ又は複数の目標構造物の中心位置とを生成するステップと、前記３Ｄ目標モデルの各々を３Ｄ球体へと変形させるステップと、変形された前記３Ｄ球体の各々に関して境界多様体を構成するステップと、前記境界多様体上でクラスタを識別し、前記１つ又は複数の目標構造物を分類するステップが行われるように構成することにより解決される。

概して、明細書中に記載された本発明の実施形態は、胸部ＣＴスキャンにおいて、小結節、血管、及び接合部のような局所構造物のタイプを分類する方法とシステムを含んでいる。この分類は肺小結節のコンピュータ支援検出（ＣＡＤ）にとって有用であり、血管や接合部により引き起こされる偽陽性（ＦＰ）を減少させるために任意の肺ＣＡＤシステムの後処理コンポーネントとして使用することができる。したがって、この分類は、ＦＰ低減の問題に注目したこのようなＣＡＤシステム又は放射線医の報告により陽性の候補が提供されることを仮定している。このようなシナリオでは、分類結果は血管や接合部により生じた偽陽性を取り除き、全体的なパフォーマンスを改善する有効な手段を提供する。

本発明の１つの実施形態による方法は様々な３Ｄトポロジカル構造の分類をはるかに単純な２Ｄデータクラスタリング分類に変換する。２Ｄデータクラスタリング分類に関しては、より包括的で柔軟な解決手段を文献から入手することができ、視覚化により適している。計算量の面での利点を別にしても、このようなアプローチは、より包括的で柔軟な分析技術のストックと、放射線医との対話の場面で有用な、より高い図解能力を持った視覚化という利点を有している。

小結節の候補が与えられると、異方性ガウシアンがデータにロバストに当てはめられる。結果として得られたガウシアンの中心及び広がりのパラメータは、当てはめられた異方性楕円体を所定サイズの異方性球体へと変形させるためにデータドメインをアフィン正規化する際に使用される。自動的な方法により、目標構造物の適切な境界面を含む３Ｄ球多様体を抽出してもよい。スケールの選択はデータ駆動式のエントロピー最小化法により行われる。この多様体から、自動モード数推定によるＥＭクラスタリング、方向統計学、及び、変形バッターチャーリヤ（Bhattacharyya）距離による階層的クラスタリングのような技術を用いて、突出構造物に対応する高輝度クラスタが分析される。推定された高輝度クラスタの個数が肺構造物のタイプを明示的に決定する：小結節（０）、付着小結節（１）、血管（２）、接合部（＞３）。本発明の１つの実施形態による方法は、方向統計学を使用して、特に多変量ラップドガウシアンモデリングを使用して、自動モード数選択を含むガウシアン当てはめ法を拡張する。

本発明の１つの実施形態による分類方法は、肺ＣＡＤの範囲を超えて、血管造影法のような様々な分野において血管構造物の有意味な情報を提供するために使用することができる。この局所的な手続きは従来技術に比べてより柔軟で効率的であり、一般的な肺ＣＡＤシステムの精度の改善に役立つであろう。さらに、モデリングの一部において選択された関心ボリューム（ＶＯＩ）表現は、ラップされていない２Ｄ境界多様体のような、ユーザ（放射線医）との対話を支援する有益な視覚化機能を有している。

胸部ＣＴイメージの選択例に関する定性的研究によれば、この分野では良好な分類結果が実証されている。本発明の１つの実施形態によるアルゴリズムは、小結節、付着小結節、血管、及び血管接合部の実例をロバストに分類することができる。

本発明は、１つの側面において、デジタル画像において肺構造物を分類する方法を提供する。この方法は、３次元（３Ｄ）デジタル画像内に１つ又は複数の目標構造物の近似的な目標構造物位置を提供するステップと、近似的な目標構造物位置に関する異方性ガウシアンモデルを当てはめ、より詳細な３Ｄ目標モデルと前記１つ又は複数の目標構造物の中心位置とを生成するステップと、３Ｄ目標モデルの各々を３Ｄ球体へと変形させるステップと、変形させた３Ｄ球体の各々に関して境界多様体を構成するステップと、前記境界多様体上でクラスタを識別し、前記１つ又は複数の目標構造物を分類するステップを有する。

本発明の別の側面によれば、デジタル画像は３次元格子上の点の領域に対応する複数の輝度値を含む。

発明の別の側面によれば、近似的な目標位置に関する異方性ガウシアンモデルの当てはめは、ガウシアンスケール空間平均シフト分析及びイェンセン-シャノン収束に基づいた自動帯域幅選択の使用と、前記目標構造物の楕円体モデルの生成とを含んでおり、前記３Ｄ楕円体モデルの中心と次元は前記ガウシアンモデルの中心と共分散に対応している。

本発明の別の側面によれば、３Ｄ目標モデルを変形させるステップは、３Ｄ楕円体をアフィン正規化するステップを含んでおり、異方性ガウシアンモデルの構造物共分散からスケーリング方向及びスケーリング係数を求める。

本発明の別の側面によれば、境界多様体を構成するステップは変形させた球体の３Ｄ表面をアンラップして２Ｄ表現を形成し、適切な境界多様体の半径を求めるステップを含む。

本発明の別の側面によれば、変形された球体の３Ｄ表面をアンラップするステップは変形された球体の表面を球座標（θ,φ）に変換するステップを含む。ただし、ここで、φ∈［-π,π］、θ∈［-π,π］である。

本発明の別の側面によれば、適切な境界多様体の半径を求めるステップは、変形された前記球体に関して異なる半径の複数の球多様体を構成するステップ、各球多様体を２Ｄ表現にアンラップするステップ、アンラップされた各球多様体上の輝度値分布を正規化するステップ、輝度値を確率として扱い、エントロピー分布を定義して、アンラップされた各球多様体の輝度エントロピーを計算するステップ、ならびに、エントロピー分布を最小化する半径を求めるステップを含む。

本発明の別の側面によれば、クラスタを識別するステップは、期待最大化法を用いて、最小記述長基準に従う前記境界多様体から突出したオブジェクトにベクトル変数
の多変量ラップドガウス分布
の混合分布
を当てはめるステップを含んでおり、前記期待最大化法の際に混合分布の要素分布の確率Ｃ_pが推定され、各次元において
は
を満たし、
は
を満たし、ｅ_k＝（0，．．．，0，1，0，．．．，0）^Tはｋ番目の成分が１でその他は０であるｋ番目のユークリッド基底ベクトルであり、混合分布の要素分布ｐの推定値
は前記期待最大化法の際に方向平均
と共分散
とに基づいて、サンプル集合
から求められ、ここで、
であり、観測値Ｘはアンラップされた２Ｄ画像Ｉ（θ,φ）から直接導かれたものであり、各サンプル値Ｉ（θ_m,φ_m）∈（-π,π］×（-π,π］の出現回数は相応する画像行列値Ｉ（θ_m,φ_m）に比例するように設定される。

本発明の別の側面によれば、上記方法は塊状の階層的クラスタリングを用いて互いに所定の距離内にあるクラスタを合併させるステップと、
に等しい多変量ラップドガウス分布の対に対して距離測度を使用するステップとをさらに有しており、ここで、μ₁およびμ₂はガウス分布の対の平均値であり、Σ₁およびΣ₂は相応する分散である。

本発明の別の側面によれば、目標構造物に関連したラップされたガウシアン要素クラスタの個数により肺構造物のクラスを決定し、その際、孤立小結節は０個のクラスタを、付着小結節は２個のクラスタを、血管は４個のクラスタを、血管接合部は６個以上のクラスタを有するものとする。

本発明の他の側面によれば、デジタル画像における肺構造物の分類のための方法ステップを実行するコンピュータにより実行可能な命令から成るプログラムを有体的に実現した、コンピュータ可読プログラム記憶装置が提供される。

ここに説明される本発明の実施形態の例は一般に胸部スキャンにおける局所構造物のタイプを分類するシステムおよび方法を含む。それに応じて、本発明は様々な変更および択一的形態を許容するが、ここでは、本発明の特定の実施形態を例として図示し、詳細に説明する。しかしながら、本発明は開示された特定の形態に限定されるものではなく、逆に、本発明の趣旨および範囲の中に入るすべての変更形態、等価形態、および択一的形態を包摂するものであることが理解されねばならない。

本明細書で使用されているように、「画像」という用語は離散的な画像要素（例えば、２Ｄ画像ならばピクセル、３Ｄ画像ならばボクセル）から成る多次元データを指す。画像は、例えば、コンピュータ断層撮影、磁気共鳴映像法、超音波、または、当業者に知られている他の任意の医用画像システムにより収集された被験者の医用画像であってよい。画像はまた、例えば、リモートセンシングシステム、電子顕微鏡などのような非医療的コンテキストから提供されるものであってもよい。画像はＲ³からＲへの関数と見なすこともできるが、本発明の方法はこのような画像に限定されるものではなく、２Ｄピクチャまたは３Ｄボリュームなどの任意の次元の画像に適用することができる。２次元または３次元画像では、画像の領域は一般には２次元または３次元の方形アレイであり、各ピクセルまたはボクセルは２つまたは３つの相互に直交する軸を基準として指定することができる。本明細書で使用される「デジタル」および「デジタル化」という用語は、適切には、デジタル収集システムを介してまたはアナログ画像からの変換を介して得られたデジタルフォーマットのまたはデジタル化されたフォーマットの画像またはボリュームを指す。

本発明の１つの実施形態による分類システムは、（１）異方性ガウシアンを当てはめるためのモジュール、（２）肺構造物の境界における２Ｄ多様体の構成、および、（３）この多様体のロバストなクラスタ分析を含む。部分（２）はエントロピー最小化に基づくデータ駆動式スケール選択を用いる。部分（３）は、自動モード数選択を用いた期待最大化法（ＥＭ）に基づくクラスタリングと、方向性データモデリングと、バッターチャーリヤ距離の変種に基づいた階層的クラスタリングを使用する。この分析における高輝度クラスタの個数が肺構造物のクラスを直接決定する。血管ツリー再構成のような他の大局的方法とは異なり、本方法は肺構造物の局所的かつ柔軟な検査を可能にする。

小結節検出に使用した場合、誤って検出されセグメント化された血管および血管の枝構造物は偽陽性（ＦＰ）の症例を意味する。本発明の１つの実施形態による分類方法はこのような非小結節構造物のすべてを棄却し、副産物として、小結節、付着小結節、血管、または血管接合部といった検討している肺構造物タイプのカテゴリーを推論する。さらに、本発明の１つの実施形態による分類解決手段は、肺ＣＡＤシステムにおける後処理フィルタとして、血管や接合部により生じたＦＰを減少させるために使用することもできる。

本発明の１つの実施形態による肺分類方法のフローチャートが図７に示されている。本発明の１つの実施形態による方法は、肺構造物の近似位置が、例えば、ＣＡＤシステム、放射線医の手動による読み取りまたはレポートなどから提供されることを仮定している。ワンクリック小結節セグメンテーションアルゴリズムを使用して、小結節、付着小結節、血管、および血管接合部を含む目標構造物を位置同定し、セグメント化してもよい。図を参照すると、ステップ７１において前もって提供された小結節の候補位置が、この半自動セグメンテーション解決手段の初期値として使用されている。ステップ７２では、異方性ガウシアンモデルが目標構造物の輝度に当てはめられ、その結果、より詳細な目標の３Ｄ楕円体モデルが得られる。これらの楕円体はステップ７３において３Ｄ球体へと変形させられる。

ステップ７４では、変形させた３Ｄ球体から境界多様体が構成される。本発明の１つの実施形態によれば、この構成は球体の３Ｄ表面を２Ｄ球座標表示へとアンラップし、それに続いて、適切な境界多様体の半径を決定することを含んでいる。境界多様体のクラスタ分析はステップ７５において行われ、続いてステップ７６において後処理が行われる。これらのステップの詳細については、以下で説明する。

ステップ７２を参照すると、本発明の１つの実施形態によるアルゴリズムは、ガウシアンスケール空間平均値シフト分析とイェンセン-シャノン収束に基づいた自動帯域幅選択とを用いて、異方性ガウシアンに基づく輝度モデルをデータにロバストに当てはめることに基づいている。これらの技術は不正確なＣＡＤまたは手動による初期化から目標物の中心の正確な推定値を提供する。３Ｄの楕円多様体は目標構造物の境界から抽出される。通常、楕円体の当てはめはトリビアルではない。しかし、このタスクはガウシアンパラメータの平均値と共分散の観点から小結節の中心の正確な推定値と楕円体形状を提供する局所構造物セグメンテーションの選択により軽減される。また、この推定のロバストネスが関心対象である血管や血管接合部／分岐のような非小結節エリアのセグメンテーションも可能にする。

数学的表現を単純化するために、元の関心ボリューム（ＶＯＩ）がステップ７３においてアフィン正規化される。このためには、セグメント化された異方性楕円体がＶＯＩの中心に位置する固定サイズの等方性球体に変換されるようにＶＯＩを変形することが必要である。アフィン正規化のパラメータ、すなわち、スケーリング方向および係数は、セグメンテーションモジュールにより推定された構造物共分散の固有値分析から直接得ることができる。

図１（ａ）−（ｃ）は、本発明の実施形態による肺構造物分類の方法を示したものである。図１（ａ）には、元々の関心ボリューム（ＶＯＩ）と、ここでは血管であるセグメント化された小結節候補とが示されており、小結節候補の小結節構造物には楕円体が当てはめられている。楕円体の当てはめは異方性ガウシアン当てはめモジュールにより得られる。図１（ｂ）には、楕円体を等方性球体に変形する元のＶＯＩのアフィン正規化が示されている。図１（ｃ）には、距離ｒ_boundにおけるセグメント化された構造物の境界多様体が示されており、この境界多様体は２Ｄ画像にアンラップされ、球座標θおよびφによりパラメータ表示されている。画像のグレイスケール値は３重線形補間を介して得られたものである。

再び図７を参照すると、ステップ７４において、検討中の肺構造物タイプのカテゴリーが、目標構造物の境界エリアから計算された適切な多様体のクラスタ分析により決定される。本発明の１つの実施形態による球多様体はアフィン正規化された３Ｄ画像から構成される。幾何学的には、この球多様体は、目標構造物の境界面から突出したオブジェクトに関する情報を含むように、この境界面を僅かに越えた球の層を表示することを目的としている。この多様体の形状は、元のＶＯＩにおいては、特に異方性ガウシアンに基づいたセグメンテーションから得られた楕円体に比例する楕円体と仮定される。したがって、アフィン正規化された表示では、この多様体は中心点（ａ_bound，ｂ_bound）と半径ｒ_boundとによって決定される等方性の球体形状に相当する。この中心点がセグメント化された楕円体の中心点と一致するのに対して、球の半径ｒ_boundは以下で説明するデータ駆動の方式により決定される。

固定されたｒ_boundを仮定すれば、境界多様体表示はデカルト座標（ｘ，ｙ，ｚ）から球座標（θ，φ）へ変換することができる。ただし、ここで、θは方位角を、φは極角を指している。その結果、アフィン正規化された楕円体の「アンラップされた」表示が２Ｄ画像行列Ｉ（θ，φ）として得られる。図１（ｃ）は、肺構造物の一例に対しての結果を示したものである。注意すべきは、一般的な慣習に反して、極角がＩｎｔｅｒｖａｌ_φ＝２π（πではなく）の区間にわたっている、すなわち、φ∈［−π，π］であり、その結果、カルテシアンボクセルが二重に出現するということである。この冗長性を導入した理由は、以下に紹介するクラスタリングでは、Ｉ（θ，φ）が両方のパラメータに関してそれぞれの区間Ｉｎｔｅｒｖａｌ_θおよびＩｎｔｅｒｖａｌ_φを越えて周期的な振る舞いを示すこと、すなわち、Ｉ（θ＋Ｉｎｔｅｒｖａｌ_θ，Ｉｎｔｅｒｖａｌ_φ）＝Ｉ（θ，φ）であることが要求されるからである。球座標の場合、Ｉｎｔｅｒｖａｌ_φ=πであれば、このことは明らかに満たされない。

本発明の１つの実施形態によれば、適切な半径ｒ_boundの決定は輝度分布のエントロピーに基づいたデータ駆動方式を用いる。この方式の動機付けのために、図２（ａ）−（ｆ）を考察する。これらの各図は、図示されているように、（θ，φ）領域におけるアンラップされた楕円体表示を異なる半径で示したものである。図２（ｇ）には、半径ｒ∈｛１，．．．，３２｝に関して、画像輝度上で計算されたそれぞれの画像ヒストグラムエントロピーが示されている。画像エントロピーは例えば
に従って画像輝度Ｉ（θ，φ）から計算することができる。

アンラップされた多様体画像は、画像輝度値分布を適切に正規化した後、２Ｄ尤度関数として扱われる。そして、確率値として解釈された正規化された輝度値から輝度エントロピーが直接計算される。半径選択は、突出構造物に起因する高輝度クラスタが相応する多様体内に最もはっきりと現れるように自動的に半径を選ぶことを必要とする。図２（ｄ）に示されているような少数のクラスタから成るこのような多様体画像は、以下の直観的な議論から、半径のより小さい画像と半径のより大きい画像よりも低いエントロピーを有するだろう。半径がより小さければ、対応する境界楕円体は目標構造物の内部を通るので、その結果、図２（ａ）−（ｂ）に示されているように尤度のより平坦な高エントロピー値が得られる。一方、半径がより大きくても、図２（ｅ）−（ｆ）に示されているように、近くにある他の「非目標」構造物の出現に起因してエントロピーは高くなる。それゆえ、適切な半径ｒ_boundがエントロピー分布Ｅ_rの局所極小値を形成する。この点から、ｒ_boundは正の差分商ΔＥ_r/Δｒが初めて出現するところに位置するように選ばれている、すなわち、

３Ｄ肺構造物の諸部分を２Ｄ画像に変換してしまえば、十分に研究された効率的で容易に視覚化可能な２Ｄ画像分析技術を適用することができる。図１（ｃ）から分かるように、境界多様体は肺構造物分類のための貴重な情報を含んでいる。実際、定められた境界から突き出しているオブジェクトの個数に等しい高輝度クラスタの個数が肺構造物のタイプを顕わにする。本発明の１つの実施形態による分類はこの観察に基づいており、領域に関する以下の仮定を念頭に置いている：
・境界多様体内に０個のクラスタがある場合には、結合した隣接構造物の不在、したがって、孤立小結節に相当するセグメント化された構造物の存在が示唆される。
・境界多様体内に１個のクラスタがある場合には、付着構造物への単一の結合が示唆される。これは、多くの場合、肺壁のようなより大きな構造物に付着した小結節から生じる。
・２個のクラスタがある場合には、２つの結合が示唆される。これは血管の場合に最もよく見られる。
・３個より多くのクラスタがある場合には、血管接合部が示唆される。

本発明の１つの実施形態によれば、高輝度クラスタの個数は図７のステップ７５において行われるクラスタリング分析を介して同定される。本発明の１つの実施形態のクラスタリング戦略は期待最大化法（ＥＭ）に基づいたガウシアンの当てはめに基づいている。標準的なＥＭガウシアンクラスタリングに加えて、本発明のクラスタリングアルゴリズムは、２Ｄ多様体画像の縁部に現れる（θ，φ）領域における連続性を反映しなければならない。特に、球の角度変数（θ，φ）をパラメータとした境界多様体表示は方向性データに対応する。方向性データを説明するために、図３の単純化された図解を考察してみよう。方向性（θ，φ）領域における適切なクラスタリングアルゴリズムは単一のクラスタを回復しなければならない。しかし、方向性モデリングの代わりに線形モデリングを使用すれば、観測可能な３つの構造物のそれぞれが独立のクラスタを形成する。さらに、本発明の１つの実施形態のクラスタリングアルゴリズムはモード数を自動的に決定できなければならない。

方向性データは超球の表面上の点として円周上に２次元で視覚化してもよい。図６（ａ）−（ｂ）には、方向性データの例が示されている。方向性データに関して生じる状況は以下の通りである。円変数θに関して、加法"ａ＋ｂ"は"（ａ＋ｂ）ｍｏｄ２π"となる。ただし、角度は区間（−π，π］内で表される。この仮定の下では、ｍｏｄ演算子も（−π，π］に写像することに注意されたい。変数
は円における対応する平均と分散を表すものとしよう。

に対する妥当な定義は、
によって表されるゼロ方向のシフトの下で不変でなければならない。円変数に対する不変量は
でなければならない。

しかし、所望の不変量が破られることは図６の例から容易に検証することができる。図６には、単純な一連の円の観察が描写されており、図６（ａ）では、Θ＝｛０．１π，０．２π，０．６π｝であり、図６（ｂ）では、
であり、これはν＝０．５πに相当する。これらの観察に関して、平均と分散の不偏最大尤度（ＭＬ）推定値は
であると計算することができる。これらは明らかにシフト不変性を犯している。図中では、平均と分散の値はそれぞれ黒点の位置と付随する円弧の長さとによって表されている。したがって、円データに関しては、平均と分散の線形的定義はゼロ方向に大きく依存しているが、このことは、相応しい操作性という点からは、不適切な挙動および要求である。

この状況に対処するため、円変数θは確率密度関数ｐ（θ）を有する確率変数であると仮定する。標準的な統計学的性質に従って、確率密度関数はｐ（θ）＞０と
を満たす。変数θは複素数ｅ^ｉθとして表示され、
により定義される円の平均方向
と円の分散
の記法を使用する。ただし、
。量ρ_θは結果長と呼ばれる。比喩的に言えば、
は期待される位相であり、ρ_θはｅ^ｉθの期待される長さである。

は散らばりの量を測るものである。容易に示されるように、平均および分散のこれらの定義は所望のシフト不変性を満たし、線形の平均および分散に対する相応しい対応物として使用することができる。

方向性データの統計的モデリングに関しては、多くのモデルが提案されている。本発明の１つの実施形態によれば、ラップドガウス分布の拡張である多変量ラップドガウス分布が使用される。線上の変数ｘのガウス分布Ｎ（ｘ）は単位半径の円の円周の周りに「巻き付け」（wrapped）ことができる。すなわち、ラップされた変数
のラップドガウス分布Ｎ_w（θ）は
である。

ベクトル変数
の多変量ラップドガウス分布は同様に
と定義することができる。ここで、ｅ_k＝（０，．．．，０，１，０，．．．，０）^Ｔはｋ番目のユークリッド基底ベクトルであり、ｋ番目の成分が１で、その他は０である。図３には、２次元多変量ラップドガウス分布の例が示されている。

方向変数の分散が適度に小さければ、
を用いてサンプル集合
から式（１）に関する正確な平均と共分散の推定値
が得られることが示されている。ここで、ｉ²＝−１であり、"ａｒｇ"は複素数の位相である。簡単のため、θのすべての次元ｆに対して２πの周期と
の変域を仮定している。

期待最大化（ＥＭ）アルゴリズムは、観測されない潜在変数にモデルが依存する確率的モデルにおけるパラメータの最大尤度推定値を求める統計的手順の１つのクラスである。ＥＭは、あたかも観測されたかのように潜在変数を含めることにより尤度の期待値を計算する予想（Ｅ）ステップと、Ｅステップにおいて求められた期待尤度を最大化することによりパラメータの最大尤度推定値を計算する最大化（Ｍ）ステップを交互に実行する。Ｍステップで求められたパラメータは新たなＥステップを開始するのに用いられ、プロセスが繰り返される。ＥＭアルゴリズムは初期推定値θ₀を繰り返し改善し、新たな推定値θ₁，．．．，θ_nを形成する。

ｙが観測変数の値から成る不完全データを表し、ｘが欠けているデータを表すものとすれば、ｘとｙとで完全なデータセットが形成される。ｐをベクトルθにより与えられるパラメータ付きの完全データの結合確率分布関数ｐ（ｙ，ｘ｜θ）とする。この関数は完全データの尤度を提供する。そこで、ベイズの公式を用いれば、観測されない変数の条件付き分布が与えられたときの期待値は
である。

この定式化は、観測不能データｐ（ｙ｜ｘ，θ）と観測可能データｐ（ｘ｜θ）が与えられたときの観測尤度に関する情報しか必要としない。θ_nからθ_n+1を導出する個々の最大化ステップは
であり、ここで、Ｅ_x［ ］はθ_nにおいて固定されたｘの条件付き分布におけるθにおけるlog ｐ（ｙ，ｘ｜θ）の条件付き期待値を表している。対数尤度ｌｏｇ ｐ（ｙ，ｘ｜θ）は真の尤度ｐ（ｙ，ｘ｜θ）の代わりに使用されることが多い。というのも、それにより式が簡単になり、しかも真の尤度と同じポイントで最大値に達するからである。言い換えれば、θ_n+1は以前のパラメータ値の下で観測変数が与えられたときの完全データの対数尤度の条件付き期待値（Ｅ）を最大化（Ｍ）する値である。一般に、最大値は対数尤度のラグランジアン関数を形成し、平均と共分散とに関して導関数を評価することにより求められる。

ＥＭクラスタリングアルゴリズムの脈絡では、正規の線形ガウシアンモデルを上に素描した多変量ラップドガウシアンモデルで置き換えることができる。特に、ＥステップとＭステップにおける元の線形バージョンがそれぞれ式（１）と、式（２）および（３）で置き換えられる。

本発明の１つの実施形態によれば、モード数選択のモデルは、ガウシアン推定の有限混合分布と、この混合分布の要素分布の個数を最小化する最小記述長（ＭＤＬ）に従うモデル選択とを使用する変形ＥＭクラスタリングアルゴリズムから求めることができる。一般に、ＥＭクラスタリングアルゴリズムへの入力は観測のサンプル集合Ｘ＝｛（θ₁，φ₁），．．．，（θ_M，φ_M）｝であり、一方で現在データは２Ｄ（画像）行列Ｉ（θ，φ）である。この非両立性を克服するため、観測ＸはＩ（θ，φ）から直接引き出され、各サンプル（θ_m，φ_m）∈（−π，π］×（−π，π］の出現回数は対応する画像行列の値Ｉ（θ_m，φ_m）に比例して設定される。

ガウシアンＥＭクラスタリングの１つの懸念は、境界多様体に突き出している構造物の真の形状のうちの１つが楕円体ガウシアンの形状に相応していない場合に生じる。このような場合には、ＥＭアルゴリズムがこの構造物にガウス分布の要素分布の集合を当てはめることが予想される。このような作用は、要素分布の個数が積分の役割を果たす分類に対して明らかに不利な影響を及ぼす。

再び図７を参照すると、本発明の１つの実施形態に従って、適当な要素分布を合併する後処理がステップ７６において行われている。特に、本発明の１つの実施形態によるこの後処理は、ＥＭで当てはめられたガウス分布のすべての要素分布から成る集合を分析し、あるスケールになるまで部分集合を単一のクラスタに合併させる第２のクラスタ分析と見なすことができる。このような状況に対する当業者に周知の１つの技術は塊状の階層的クラスタリングである。階層的クラスタリングでは、クラスタ空間はその要素の距離によって表現される。今のケースでは、要素は多変量ラップドガウシアン関数であり、幾何学形状に対して統計的記述子が使用される。ガウス分布に相応しい（さらに解析的に計算可能な）統計的距離測度はバッターチャーリヤ距離
である。しかし、Ｄ_Bhattはラップドガウシアンの方向特性を考慮しない。したがって、本発明の１つの実施形態によれば、Ｄ_Bhattの変種である「ラップドバッターチャーリヤ距離」
が提案される。最終的に、ラップドガウシアンの要素クラスタの個数が肺構造物のクラスを決定する：０であれば、孤立小結節であり、２×１＝２であれば、付着小結節であり、２×２＝４であれば、血管であり、＞２×３＝６であれば、血管接合部である。係数２は上で論じた極座標φにおける二重区間に起因するものである。

本発明の１つの実施形態による方法の限界は、（θ，φ）領域内においてスケールが位置依存することである。本発明の１つの実施形態による１つ代替形態では、方向性データをフォン・ミーゼス-フィッシャー分布でモデリングすることでこの問題が回避される。ｄ次元のランダム単位ベクトルに対するフォン・ミーゼス-フィッシャー分布は次式の形をとる。

f(x|μ,κ)=c_d(κ)exp(κμ^Tx),
ここで、‖μ‖＝１、κ≧０、ｄ≧２である。正規化定数ｃ_d（κ）は
により与えられ、ここで、Ｉ_r（ ）はｒ次の変形第１種ベッセル関数を表している。分布は平均方向μと濃度パラメータκとによりパラメトライズされている。なお、濃度パラメータκは、ｆ（ｘ｜μ，κ）に従って引かれた単位ベクトルが平均方向μの周りにどの程度集中しているかを特徴付けるものである。しかし、フォン・ミーゼス-フィッシャー分布に関するパラメータ推定のためには、一般に解析解の存在しないベッセル関数の比の陰方程式を解く必要がある。

本発明の１つの実施形態によれば、提案された肺構造物分類に関する定性的な実験が行われた。図４および５は胸部ＣＴ画像の分類の図を示したものであり、「小結節」、「付着小結節」、「血管」、「血管接合部」の各クラスについて２つの例が示されている。

図４（ａ）−（ｄ）および５（ａ）−（ｄ）は、胸部ＣＴスキャンに関して、本発明の１つの実施形態の肺構造物分類方法の図解例を示したものである。各行は１つの例のセグメンテーションと検証とに対応している。図４の最初の２行は小結節オブジェクトに関するものであり、最後の２行は肺壁に付着した小結節に関するものであり、一方、図５の行１および２は血管構造物を、行３および４は血管接合部を示している。列（ａ）はＣＴ ＶＯＩを３つの直交する断面で示したものである。セグメンテーション結果は楕円により表されている。列（ｂ）は、３Ｄ楕円体が球へと変形するような元のＶＯＩのアフィン正規化を描写している。列（ｃ）は（θ,φ）領域にアンラップされた構成された境界多様体を示している。しかし、付加的な輝度閾が導入されていることに注意されたい。このステップはガウシアンEMクラスタリングを混乱させかねない低輝度構造物を除去する迅速かつ簡単な手段として使用される。列（ｄ）の図はEMに基づいたアルゴリズムによるガウシアン混合分布モデルの当てはめの結果を示している。破線の楕円はEMに基づいてクラスタリングされたがウシアン要素分布に対応しており、実線の楕円は後処理の後のクラスタを描写したものである。

列（ａ）に示されているように、３Ｄセグメンテーションは、図４に示されているようにすべての孤立小結節及び付着小結節をセグメント化することも、図５に示されているように偽陽性の血管及び血管接合部をセグメント化することもできる。列（ｄ）では、境界多様体画像がセクション２．２．１で説明したようにサンプルデータ集合Ｘに変換されている。さらに、列（ｄ）はＥＭに基づいたラップドガウシアンクラスタリングの結果を示している、すなわち、ｋ個の要素分布の平均と共分散が破線の楕円により表されている。特に、図４の行３及び４と図５の行３及び４の（θ,φ）領域の縁部における連続性に注意されたい。視覚化目的のために、階層的クラスタリング後処理の実例が含まれている。この後処理からのクラスタはk₂個の実線の楕円により表されており、この楕円の中心点と広がりは、後処理された１つのクラスタの中のすべてのラップドガウシアンの平均から計算された平均と共分散とに対応している。クラスタの濃度が小さければ、この実例は、例えば図５の行２におけるような、退化した楕円に導くことがありうる。要素分布の個数k₂から構造物のクラスを推論することにより、提示された分類が８つのすべての例に対して正しい答えを与えているかを検証することができる。他のケースでも同様の結果が得られた。

分類の限界を指摘することは価値のあることである。というのも、ある状況においては、誤った分類に導かれかねないからである。３Ｄ球多様体の極（θ＝０とφ＝πに相当する）における構造物はアンラップ後の２Ｄ画像のθ方向において不釣合いに大きくなる。この状況は、北極地域と南極地域が３Ｄ球体の地球儀におけるよりも２Ｄメルカトル投影された世界地図において比較的大きな地域を占めるような地図作製法における現象と比較することができる。上記の例では、この挙動は
における高輝度構造物がθの（−π，π］の全区間を超えて延びている図５の行４に見られる。その結果、アンラップされた多様体におけるスケール関係、特にこれらの極領域に対するスケール関係から結論を導く際には、注意が必要である実際、これはラップドガウシアンモデリングの、特にアンラッピングの、欠点である。この点、フォン・ミーゼス-フィッシャーモデリングではアンラッピングが仮定されないので、この現象が回避されることに注意されたい。

本発明は様々な形態のハードウェア、ソフトウェア、ファームウェア、特殊用途向けプロセス、又はこれらの組合せにより実施可能であることが理解されなければならない。１つの実施形態では、本発明をコンピュータ可読プログラム記憶装置上に有体的に実現されたアプリケーションプログラムとしてソフトウェアで実施することもできる。このアプリケーションプログラムは任意の相応しいアーキテクチャを持つマシーンにアップロードし、実行することができる。

図８は、本発明の１つの実施形態による分類方法を実現するコンピュータシステムの一例のブロック図である。図８を参照すると、本発明を実施するコンピュータシステム８１は、とりわけ、中央処理ユニット（ＣＰＵ）８２、メモリ８３、及び入出力（Ｉ／Ｏ）インタフェース８４を有する。コンピュータシステム８１は一般にＩ／Ｏインタフェース８４を介してディスプレイ８５や、マウス及びキーボードのような様々な入力デバイス８６と結合されている。補助回路としては、キャッシュ、電源、クロック回路、及び通信バスのような回路を含んでいてよい。メモリ８３には、ランダムアクセスメモリ（ＲＡＭ）、読出し専用メモリ（ＲＯＭ）、ディスクドライブ、テープドライブなど、またはこれらの組合せが含まれうる。本発明は、メモリ８３に格納され、ＣＰＵ８２により実行されることにより、信号源８８からの信号を処理するルーチン８７として実施することができる。コンピュータシステム８１はそれ自体としては汎用コンピュータシステムであり、本発明のルーチン８７を実行する際には専用コンピュータシステムとなる。

また、コンピュータシステム８１はオペレーティングシステムとマイクロ命令コードも含んでいる。本明細書に記載された様々なプロセスや機能はマイクロ命令コードの一部であってもよいし、オペレーティングシステムを介して実行されるアプリケーションプログラム（またはこれらの組合せ）の一部であってもよい。加えて、他の様々な周辺機器を、付加的なデータ記憶装置や印刷機のようなコンピュータプラットフォームに接続してもよい。

さらに、構成要素であるシステムコンポーネントと添付図面に描かれた方法ステップのうちのいくつかはソフトウェアとして実施することができるので、システムコンポーネント（又はプロセスステップ）間の実際の接続は本発明のプログラムの仕方に応じて異なりうる。本明細書に開示された本発明の教示に従えば、当業者は本発明のこれらの実施形態又は構成ならびにこれらと類似した実施形態及び構成を考えることができる。

有利な実施形態を参照して本発明を詳細に説明してきたが、当業者には、添付された請求項において示されている本発明の趣旨及び範囲から逸脱することなく、本発明に対して様々な変更及び置き換えを為しうることが理解されるであろう。

（ａ）−（ｃ）は、本発明の１つの実施形態に従って肺構造物分類の方法を示す。 （ａ）−（ｇ）は、本発明の１つの実施形態に従って、半径の異なるアンラップされた楕円体の効果と相応する画像輝度ヒストグラムのエントロピーとを示す。 本発明の１つの実施形態に従って、方向性データを用いたクラスタリングを示す。 （ａ）−（ｄ）は、胸部ＣＴスキャンに関して、本発明の１つの実施形態の肺構造物分類方法の図解例を示す。 （ａ）−（ｄ）は、胸部ＣＴスキャンに関して、本発明の１つの実施形態の肺構造物分類方法の図解例を示す。 （ａ）−（ｂ）は、本発明の１つの実施形態に従って、方向性データの例を示す。 本発明の１つの実施形態による分類方法のフローチャートである。 本発明の１つの実施形態による分類方法を実現するコンピュータシステムの一例のブロック図である。

Claims (24)

1. デジタル画像において肺構造物を分類する方法において、３次元（３Ｄ）デジタル画像内の１つ又は複数の目標構造物の近似的な目標構造物位置を提供するステップと、前記の近似的な目標構造物位置に関して異方性ガウシアンモデルを当てはめ、より詳細な３Ｄ目標モデルと前記１つ又は複数の目標構造物の中心位置とを生成するステップと、前記３Ｄ目標モデルの各々を３Ｄ球体へと変形させるステップと、変形された前記３Ｄ球体の各々に関して境界多様体を構成するステップと、前記境界多様体上でクラスタを識別し、前記１つ又は複数の目標構造物を分類するステップを有することを特徴とする、デジタル画像において肺構造物を分類する方法。
2. 前記デジタル画像は３次元格子上の点の領域に対応する複数の輝度値を含む、請求項１記載の方法。
3. 近似的な目標位置に関して異方性ガウシアンモデルの当てはめる前記ステップは、ガウシアンスケール空間平均シフト分析とイェンセン-シャノン収束に基づいた自動帯域幅選択とを使用するステップと、前記目標構造物の楕円体モデルを生成し、該３Ｄ楕円体モデルの中心と次元を前記ガウシアンモデルの中心と共分散に対応させるステップを含む、請求項１記載の方法。
4. ３Ｄ目標モデルを変形させる前記ステップは、前記３Ｄ楕円体をアフィン正規化し、前記異方性ガウシアンモデルの構造物共分散からスケーリング方向及びスケーリング係数を求めるステップを含む、請求項３記載の方法。
5. 境界多様体を構成する前記ステップは、変形された前記球体の３Ｄ表面をアンラップして２Ｄ表現を形成し、適切な境界多様体の半径を求めるステップを含む、請求項１記載の方法。
6. 変形された前記球体の３Ｄ表面をアンラップする前記ステップは変形された前記球体の表面を球座標(θ,φ)に変換するステップを含み、ここで、φ∈[-π,π]かつθ∈[-π,π]である、請求項５記載の方法。
7. 適切な境界多様体の半径を求める前記ステップは、変形された前記球体に関して異なる半径の複数の球多様体を構成するステップと、各球多様体を２Ｄ表現へアンラップするステップと、アンラップされた前記各球多様体上の輝度値分布を正規化するステップと、輝度値を確率として扱い、エントロピー分布を定義して、アンラップされた前記各球多様体の輝度エントロピーを計算するステップと、前記エントロピー分布を最小化する半径を求めるステップを含む、請求項５記載の方法。
8. クラスタを識別する前記ステップは、期待最大化法を用いて、最小記述長基準に従う前記境界多様体から突出したオブジェクトにベクトル変数
   のラップされた多変量ガウス分布
   の混合分布
   を当てはめるステップを含んでおり、ただし、混合分布の要素分布の確率Ｃ_pは前記期待最大化法の際に推定されたものであり、各次元において
   は
   を満たし、
   は
   を満たし、ｅ_k＝（0，．．．，0，1，0，．．．，0）^Tはｋ番目の成分が１でその他は０であるｋ番目のユークリッド基底ベクトルであり、混合分布の要素分布ｐの推定値
   は前記期待最大化法の際に方向平均
   と共分散
   とに基づいて、サンプル集合
   から求められ、ここで、
   であり、観測値Ｘはアンラップされた２Ｄ画像Ｉ（θ,φ）から直接導かれたものであり、各サンプル値Ｉ（θ_m,φ_m）∈（-π,π］×（-π,π］の出現回数は相応する画像行列値Ｉ（θ_m,φ_m）に比例するように設定される、請求項１記載の方法。
9. 塊状の階層的クラスタリングを用いて互いに所定の距離内にあるクラスタを合併させるステップと、
   に等しいラップされた多変量ガウス分布の対に対して距離測度を使用するステップとをさらに有しており、ここで、μ₁およびμ₂はガウス分布の対の平均値であり、Σ₁およびΣ₂は相応する分散である、請求項１記載の方法。
10. 目標構造物に関連したラップされたガウシアン要素クラスタの個数により肺構造物のクラスを決定し、その際、孤立小結節は０個のクラスタを、付着小結節は２個のクラスタを、血管は４個のクラスタを、血管接合部は６個以上のクラスタを有するものとする、請求項１記載の方法。
11. デジタル画像において肺構造物を分類する方法において、３次元格子上の点の領域に対応する複数の輝度を有する３次元（３Ｄ）デジタル画像内に、該画像内の目標構造物を表す１つまたは複数の３Ｄ球体の目標位置を提供するステップと、前記３Ｄ球体に関して半径の異なる複数の球多様体を構成するステップと、輝度値を確率値として扱い、エントロピー分布を定めて、前記各球多様体の輝度分布を計算するステップと、前記エントロピー分布を最小化する半径を求め、この半径最小化により境界多様体を定めるステップと、境界多様体の表面を２Ｄ球座標（θ，φ）表示へとアンラップするステップと、ただし、φ∈［−π，π］かつθ∈［−π，π］であり、期待最大化法を用いてベクトル変数
    のラップされた多変量ガウス分布Ｎ_w（Θ）の混合分布
    を前記境界多様体から突出している目標構造物のクラスタに当てはめるステップを有しており、ただし、混合分布の要素分布の確率Ｃ_pは前記期待最大化法の際に推定されたものであり、
    であり、肺構造物は突出しているクラスタの個数により分類される、ことを特徴とするデジタル画像において肺構造物を分類する方法。
12. 前記複数の球多様体の各々における輝度分布を正規化するステップをさらに有する、請求項１１記載の方法。
13. １つまたは複数の３Ｄ球体の目標位置を提供する前記ステップは、前記３次元（３Ｄ）デジタル画像内における前記１つまたは複数の目標構造物の近似的な目標構造物位置を提供するステップと、前記近似的目標位置に関して異方性ガウシアンモデルを当てはめ、より詳細な３Ｄ楕円体モデルと該３Ｄ楕円体モデルの中心位置とを生成するステップと、前記楕円体モデルを３Ｄ球体へとアフィン正規化するステップを含む、請求項１１記載の方法。
14. 各次元において
    は
    を満たし、Ｎ_w（Θ）は
    を満たし、ｅ_k＝（0，．．．，0，1，0，．．．，0）^Tはｋ番目の成分が１でその他は０であるｋ番目のユークリッド基底ベクトルであり、混合分布の要素分布ｐの推定値
    は前記期待最大化法の際に方向平均
    と共分散
    とに基づいて、サンプル集合
    から求められ、ここで、
    であり、観測値Ｘはアンラップされた２Ｄ画像Ｉ（θ,φ）から直接導かれたものであり、各サンプル値Ｉ（θ_m,φ_m）∈（-π,π］×（-π,π］の出現回数は相応する画像行列値Ｉ（θ_m,φ_m）に比例するように設定される、請求項１１記載の方法。
15. デジタル画像内の肺構造物を分類する方法ステップを実行するコンピュータにより実行可能な命令から成るプログラムを有体的に実現した、コンピュータ可読プログラム記憶装置において、前記方法が、デジタル化された３次元（３Ｄ）画像内の１つ又は複数の目標構造物の近似的な目標構造物位置を提供するステップと、前記の近似的な目標構造物位置に関する異方性ガウシアンモデルを当てはめ、より詳細な３Ｄ目標モデルと前記１つ又は複数の目標構造物の中心位置とを生成するステップと、前記３Ｄ目標モデルの各々を３Ｄ球体へと変形させるステップと、変形させた前記３Ｄ球体の各々に関して境界多様体を構成するステップと、前記境界多様体上でクラスタを識別し、前記１つ又は複数の目標構造物を分類するステップを有することを特徴とする、コンピュータ可読プログラム記憶装置。
16. 前記デジタル画像は３次元格子上の点の領域に対応する複数の輝度値を含む、請求項１5記載のコンピュータ可読プログラム記憶装置。
17. 近似的な目標位置に関する異方性ガウシアンモデルの当てはめは、ガウシアンスケール空間平均シフト分析及びイェンセン-シャノン収束に基づいた自動帯域幅選択の使用と、前記目標構造物の楕円体モデルの生成とを含んでおり、前記３Ｄ楕円体モデルの中心と次元は前記ガウシアンモデルの中心と共分散に対応している、請求項１５記載のコンピュータ可読プログラム記憶装置。
18. ３Ｄ目標モデルを変形させる前記ステップは前記３Ｄ楕円体をアフィン正規化するステップを含んでおり、前記異方性ガウシアンモデルの構造物共分散からスケーリング方向及びスケーリング係数が求められる、請求項１７記載のコンピュータ可読プログラム記憶装置。
19. 境界多様体を構成する前記ステップは変形させた球体の３Ｄ表面を２Ｄ表現へとアンラップし、適切な境界多様体の半径を求めるステップを含む、請求項１５記載のコンピュータ可読プログラム記憶装置。
20. 変形された前記球体の３Ｄ表面を展開する前記ステップは変形された前記球体の表面を球座標（θ，φ）に変換するステップを含み、ここで、φ∈［−π，π］かつθ∈［−π，π］である、請求項１９記載のコンピュータ可読プログラム記憶装置。
21. 適切な境界多様体の半径を求める前記ステップは、変形された前記球体に関して異なる半径の複数の球多様体を構成するステップと、各球多様体を２Ｄ表現へとアンラップするステップと、アンラップされた前記各球多様体上の輝度値分布を正規化するステップと、輝度値を確率として扱い、エントロピー分布を定義して、アンラップされた前記各球多様体の輝度エントロピーを計算するステップと、前記エントロピー分布を最小化する半径を求めるステップを含む、請求項１９記載のコンピュータ可読プログラム記憶装置。
22. クラスタを識別する前記ステップは、期待最大化法を用いて、最小記述長基準に従う前記境界多様体から突出したオブジェクトにベクトル変数
    のラップされた多変量ガウス分布
    の混合分布
    を当てはめるステップを含んでおり、前記期待最大化法の際に混合分布の要素分布の確率Ｃ_pが推定され、各次元において
    は
    を満たし、
    は
    を満たし、ｅ_k＝（0，．．．，0，1，0，．．．，0）^Tはｋ番目の成分が１でその他は０であるｋ番目のユークリッド基底ベクトルであり、混合分布の要素分布ｐの推定値
    は前記期待最大化法の際に方向平均
    と共分散
    とに基づいて、サンプル集合
    から求められ、ここで、
    であり、観測値Ｘはアンラップされた２Ｄ画像Ｉ（θ,φ）から直接導かれたものであり、各サンプル値Ｉ（θ_m,φ_m）∈（-π,π］×（-π,π］の出現回数は相応する画像行列値Ｉ（θ_m,φ_m）に比例するように設定される、請求項１５記載のコンピュータ可読プログラム記憶装置。
23. 前記方法はさらに、塊状の階層的クラスタリングを用いて互いに所定の距離内にあるクラスタを合併させるステップと、
    に等しい多変量ラップドガウス分布の対に対して距離測度を使用するステップとをさらに有しており、ここで、μ₁およびμ₂はガウス分布の対の平均値であり、Σ₁およびΣ₂は相応する分散である、請求項１５記載のコンピュータ可読プログラム記憶装置。
24. 目標構造物に関連したラップドガウシアンの要素クラスタの個数により肺構造物のクラスを決定し、その際、孤立小結節は０個のクラスタを、付着小結節は２個のクラスタを、血管は４個のクラスタを、血管接合部は６個以上のクラスタを有するものとする、請求項１５記載のコンピュータ可読プログラム記憶装置。