(57)【特許請求の範囲】
【請求項１】
ＣＴ値に対応する閾値に基づいて、３次元画像データから相対的に輝度値が高く、抽出対象の構造を含む第１の領域と、相対的に輝度値が低く、前記３次元画像データにおける空気領域である第２の領域とを抽出する第１の抽出部と、
前記第１の領域をボクセルの特徴量の類似度及び位置関係に基づいて複数の小領域に分割し、前記第２の領域に隣接し、前記第２の領域との類似度が相対的に高い小領域を前記複数の小領域から除外することで、前記複数の小領域から前記構造に含まれる小領域を抽出して、前記構造の領域を抽出する第２の抽出部と、
を備える画像処理装置。
【請求項２】
前記第２の抽出部は、前記第２の領域と所定の位置関係にある小領域のうち、前記第２の領域との特徴量分布の類似度が所定の閾値以下、又は、前記特徴量分布の距離が所定の閾値以上である小領域を前記構造に含まれる小領域として抽出する請求項１に記載の画像処理装置。
【請求項３】
前記第２の抽出部は、前記抽出対象の構造に該当するかを特定するための特定情報に基づいて、前記複数の小領域から前記構造に含まれる小領域を抽出する請求項１又は２に記載の画像処理装置。
【請求項４】
前記第２の抽出部は、前記第１の領域に含まれるボクセルを輝度値及び３次元座標を用いてクラスタリングすることにより、前記第１の領域を複数の小領域に分割する請求項１〜３のいずれか一項に記載の画像処理装置。
【請求項５】
前記特徴量分布がボクセルの輝度値の分布である請求項２に記載の画像処理装置。
【請求項６】
前記第１の抽出部は、前記３次元画像データにおいてボクセルの輝度値が所定の閾値以上である領域を前記第１の領域として抽出する請求項１〜５のいずれか一項に記載の画像処理装置。
【請求項７】
前記３次元画像データがＣＴ画像データであり、前記抽出対象の構造が腫瘍及び血管を含む解剖学的構造である請求項１〜６のいずれか一項に記載の画像処理装置。
【請求項８】
前記腫瘍は、充実性陰影及びすりガラス状陰影のうち、少なくとも一方の陰影として観察される請求項７に記載の画像処理装置。
【請求項９】
前記構造の領域に含まれる小領域を識別する識別部をさらに備える請求項１〜８のいずれか一項に記載の画像処理装置。
【請求項１０】
前記特定情報が、操作者によって指定された情報、又は、ＣＡＤによって指定された情報である請求項３に記載の画像処理装置。
【請求項１１】
３次元画像データから抽出対象の領域を抽出する装置によって実行される画像処理方法であって、
ＣＴ値に対応する閾値に基づいて、前記３次元画像データから相対的に輝度値が高く、抽出対象の構造を含む第１の領域と、相対的に輝度値が低く、前記３次元画像データにおける空気領域である第２の領域とを抽出する第１の抽出ステップと、
前記第１の領域をボクセルの特徴量の類似度及び位置関係に基づいて複数の小領域に分割し、前記第２の領域に隣接し、前記第２の領域との類似度が相対的に高い小領域を前記複数の小領域から除外することで、前記複数の小領域から前記構造に含まれる小領域を抽出して、前記構造の領域を抽出する第２の抽出ステップと、
を含む画像処理方法。
【請求項１２】
３次元画像データを収集する画像データ収集部と、
ＣＴ値に対応する閾値に基づいて、前記３次元画像データから相対的に輝度値が高く、抽出対象の構造を含む第１の領域と、相対的に輝度値が低く、前記３次元画像データにおける空気領域である第２の領域とを抽出する第１の抽出部と、
前記第１の領域をボクセルの特徴量の類似度及び位置関係に基づいて複数の小領域に分割し、前記第２の領域に隣接し、前記第２の領域との類似度が相対的に高い小領域を前記複数の小領域から除外することで、前記複数の小領域から前記構造に含まれる小領域を抽出して、前記構造の領域を抽出する第２の抽出部と、
を備える医用画像診断装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、画像処理装置、画像処理方法及び医用画像診断装置に関する。
【背景技術】
【０００２】
従来、Ｘ線ＣＴ（Computed  Tomography）装置などの医用画像診断装置によって収集された３次元の画像データに含まれる関心領域を自動で抽出する技術が知られており、例えば、腫瘍の鑑別診断に利用されている。一例を挙げると、Ｘ線ＣＴ装置によって異なる時期に収集された複数の画像データに含まれる結節領域をそれぞれ抽出して、抽出した結節領域の体積や直径の変化を調べることにより、結節が悪性であるか又は良性であるかを鑑別する。
【０００３】
このような鑑別診断は、例えば、肺がん等の疑いを持つ患者に対して行われ、近年では、結節領域を抽出するだけではなく、結節領域をすりガラス状の半透明の画像領域であるＧＧＯ  (Ground  Glass  Opacity)  領域や、ＧＧＯよりも高輝度な画像領域であるソリッド(Solid:充実性)領域などに分類して、それぞれの体積や直径の変化により鑑別診断が行われている。かかる鑑別診断の精度を向上させるためには、関心領域を正確に抽出することが重要である。
【０００４】
現在、３次元の画像データから関心領域を抽出する方法としては、例えば、ユーザ又はＣＡＤ（Computer-Aided-Diagnosis）システムによって指定された腫瘍候補点の周囲に３次元のＲＯＩ（Region  Of  Interest）であるＶＯＩ（Volume  Of  Interest）を設定して、設定したＶＯＩに含まれるボクセルを輝度の閾値処理によって区分けすることで腫瘍候補領域を抽出する。そして、抽出した腫瘍候補領域から血管などの領域を除去して、Ｓｏｌｉｄ領域及びＧＧＯ領域に組成分類する方法が知られている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２００９−２８１６１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
本発明が解決しようとする課題は、３次元画像データからの領域抽出の精度を向上させることができる画像処理装置、画像処理方法及び医用画像診断装置を提供することである。
【課題を解決するための手段】
【０００７】
実施形態の画像処理装置は、第１の抽出部と、第２の抽出部とを備える。第１の抽出部は、ＣＴ値に対応する閾値に基づいて、３次元画像データから相対的に輝度値が高く、抽出対象の構造を含む第１の領域と、相対的に輝度値が低く、前記３次元画像データにおける空気領域である第２の領域とを抽出する。第２の抽出部は、前記第１の領域をボクセルの特徴量の類似度及び位置関係に基づいて複数の小領域に分割し、前記第２の領域に隣接し、前記第２の領域との類似度が相対的に高い小領域を前記複数の小領域から除外することで、前記複数の小領域から前記構造に含まれる小領域を抽出して、前記構造の領域を抽出する。
【図面の簡単な説明】
【０００８】
【図１】第１の実施形態に係る医用画像処理システムの一例を示す図。
【図２】ボリュームデータからの領域抽出の一例を説明するための図。
【図３】従来技術に係る課題を説明するための図。
【図４】第１の実施形態に係る画像処理装置の一例を示す図。
【図５】第１の実施形態に係る画像処理装置による処理の手順を示すフローチャート。
【図６Ａ】第１の領域の抽出の一例を示す図。
【図６Ｂ】閾値の決定の一例を説明するための図。
【図７】第１の実施形態に係る第２の抽出部による処理の概念を説明するための模式図。
【図８Ａ】第１の実施形態に係る第２の抽出部による領域分割の一例を示す図。
【図８Ｂ】第１の実施形態に係る第２の抽出部によるクラスタリングの一例を説明するための模式図。
【図９】第１の実施形態に係る第２の抽出部による抽出処理の一例を示す図。
【図１０】第２の実施形態に係る画像処理装置の一例を示す図。
【図１１】第２の実施形態に係る画像処理装置による処理の手順を示すフローチャート。
【図１２】第２の実施形態に係る識別部による処理の一例を説明するための図。
【図１３】実施形態に係る画像処理装置のハードウェア構成を示す図。
【発明を実施するための形態】
【０００９】
以下、添付図面を参照して、実施形態に係る画像処理装置、画像処理方法及び医用画像診断装置を説明する。以下の実施形態では、同一の参照符号を付した部分は同様の動作を行うものとして、重複する説明を適宜省略する。
【００１０】
（第１の実施形態）
まず、第１の実施形態に係る画像処理装置を含む医用情報処理システムについて説明する。図１は、第１の実施形態に係る医用情報処理システム１の一例を示す図である。図１に示すように、第１の実施形態に係る医用情報処理システム１は、画像処理装置１００と、医用画像診断装置２００と、画像保管装置３００とを有する。図１に例示する各装置は、例えば、病院内に設置された院内ＬＡＮ（Local  Area  Network）により、直接的、又は間接的に相互に通信可能な状態となっている。例えば、医用情報処理システム１にＰＡＣＳ（Picture  Archiving  and  Communication  System）が導入され、各装置は、ＤＩＣＯＭ（Digital  Imaging  and  Communications  in  Medicine）規格に則って、３次元の画像データ（ボリュームデータ）等を相互に送受信する。
【００１１】
医用画像診断装置２００は、Ｘ線診断装置、Ｘ線ＣＴ（Computed  Tomography）装置、ＭＲＩ（Magnetic  Resonance  Imaging）装置、超音波診断装置、ＳＰＥＣＴ（Single  Photon  Emission  Computed  Tomography）装置、ＰＥＴ（Positron  Emission  computed  Tomography）装置、ＳＰＥＣＴ装置とＸ線ＣＴ装置とが一体化されたＳＰＥＣＴ−ＣＴ装置、ＰＥＴ装置とＸ線ＣＴ装置とが一体化されたＰＥＴ−ＣＴ装置、又はこれらの装置群等である。また、第１の実施形態に係る医用画像診断装置２００は、ボリュームデータを生成可能である。
【００１２】
具体的には、医用画像診断装置２００は、被検体を撮影することによりボリュームデータを生成する。例えば、医用画像診断装置２００は、被検体を撮影することにより投影データやＭＲ信号等のデータを収集し、収集したデータから、被検体の体軸方向に沿った複数のアキシャル面の医用画像データを再構成することで、ボリュームデータを生成する。例えば、医用画像診断装置２００は、５００枚のアキシャル面の医用画像データを再構成する。この５００枚のアキシャル面の医用画像データ群が、ボリュームデータである。なお、医用画像診断装置２００により撮影された被検体の投影データやＭＲ信号等自体をボリュームデータとしても良い。
【００１３】
また、医用画像診断装置２００は、生成したボリュームデータを画像保管装置３００に送信する。なお、医用画像診断装置２００は、ボリュームデータを画像保管装置３００に送信する際に、付帯情報として、例えば、患者を識別する患者ＩＤ、検査を識別する検査ＩＤ、医用画像診断装置２００を識別する装置ＩＤ、医用画像診断装置２００による１回の撮影を識別するシリーズＩＤ等を送信する。
【００１４】
画像保管装置３００は、医用画像を保管するデータベースである。具体的には、画像保管装置３００は、医用画像診断装置２００から送信されたボリュームデータを記憶部に格納し、これを保管する。なお、本実施形態は、大容量の画像を保管可能な画像処理装置１００を用いることで、図１に例示する画像処理装置１００と画像保管装置３００とが統合される場合であっても良い。すなわち、本実施形態は、画像処理装置１００そのものにボリュームデータを記憶させる場合であっても良い。
【００１５】
なお、第１の実施形態において、画像保管装置３００に保管されたボリュームデータは、患者ＩＤ、検査ＩＤ、装置ＩＤ、シリーズＩＤ等と対応付けて保管される。このため、画像処理装置１００は、患者ＩＤ、検査ＩＤ、装置ＩＤ、シリーズＩＤ等を用いた検索を行なうことで、必要なボリュームデータを画像保管装置３００から取得する。
【００１６】
画像処理装置１００は、医用画像診断装置２００又は画像保管装置３００からボリュームデータを受信し、受信したボリュームデータに含まれる所定の領域を抽出する。ここで、本実施形態に係る画像処理装置１００は、従来技術と比較して、ボリュームデータからの領域抽出の精度を向上させる。具体的には、画像処理装置１００は、ボリュームデータから抽出対象の構造を含む領域を抽出する際の精度を向上させる。
【００１７】
ここで、ボリュームデータからの領域抽出の一例と従来技術における課題について、図２及び３を用いて説明する。図２は、ボリュームデータからの領域抽出の一例を説明するための図である。なお、図２においては、Ｘ線ＣＴ装置によって収集された胸部のＣＴ画像データから結節領域を抽出して、組成（Ｓｏｌｉｄ領域及びＧＧＯ領域）ごとの体積を推定する場合についての模式図を示す。
【００１８】
ボリュームデータからの領域抽出においては、例えば、図２に示すように、まず、「入力」として胸部のＣＴボリュームデータを受け付ける。ここで、抽出対象である結節の中心点（図中の×印）や結節が含まれる肺野領域が指定される。このとき、中心点や肺野領域の指定は、操作者によって手動で行われたり、ＣＡＤ（Computer-Aided-Diagnosis）システムによって自動で行われたりする。このように抽出対象が指定されると、抽出対象を含む前景領域が抽出される。
【００１９】
例えば、図２の「前景領域」に示すように、指定された肺野領域の中から抽出対象である結節を含む前景領域が抽出される。ここで、前景領域の抽出では、ボクセルにおける輝度値に基づいて前景領域が抽出されるため、結節と同様に輝度値が高い血管が前景領域として抽出される。そこで、次に、抽出された前景領域に含まれる血管領域が除去され（図中の「血管除去」）、結節の領域を示す境界が修正されることで（図中の「境界修正」）、ボリュームデータから結節の領域が抽出される。その後、抽出された結節の領域におけるＳｏｌｉｄ領域及びＧＧＯ領域が分類され（図中の「組成分類」）、「出力」として、それぞれの領域が示された画像とともに推定された体積値が表示される。
【００２０】
上述したように、ボリュームデータからの抽出対象の抽出には、まず、抽出対象を含む前景領域が抽出され、抽出された前景領域から対象領域が抽出される。しかしながら、従来技術においては、上述した前景領域の抽出の精度に一定の限界があるため、最終的な出力である抽出対象の組成（構造）ごとの領域表示や、体積推定の精度においても一定の限界があった。
【００２１】
例えば、ＣＴボリュームデータから結節を抽出する場合、ボクセルにおけるＣＴ値（輝度値）を閾値と比較することで前景領域が抽出される。ＣＴ値は、Ｘ線の吸収率を示し，「空気」を−１０００ＨＵ、「水」を０ＨＵと定義したＨＵ（Hounsfield  Unit）を単位として利用する。ＣＴ値は、一般に空気（−１０００）を基準として、−１０００〜１０００の値の範囲をとり、体内にある種々の解剖学的構造のＣＴ値は−８００〜３００の範囲に分布する。例えば、ＧＧＯ領域は−８００〜−５００、血管やＳｏｌｉｄ領域は−５００〜０、胸壁や骨は０〜５００に分布することが多い。従って、輝度値の閾値処理によってこれらの解剖学的構造を抽出する場合には、閾値はＧＧＯ領域が抽出可能な−８００前後に設定することが望ましい。
【００２２】
しかしながら、ＣＴボリュームデータはスライス厚や解像度などのパラメータ、再構成関数等によるボケの影響で、解剖学的構造の境界付近の輝度値が本来より高くなることがしばしばある。このようなデータに対して輝度値の閾値処理を行った場合、空気よりわずかに輝度値の高い境界付近の領域が過検出されてしまう。特に腫瘍や血管が複雑に入り組んでいる場合、近接する構造同士が過検出領域を介して結合し、後段で除去処理を行ったとしても正しく分離することが難しいという問題がある。例えば、図３に示すように、結節と血管とを含むＣＴボリュームデータに対して輝度値の閾値処理を行うと（図中の「輝度値による前景抽出」）、結節や血管の周囲の領域が過検出されてしまい（図中の「過検出領域」）、その後の「分類」で精度の高い分類を行うことが困難となる。このような問題に対しては、対策として閾値を引き上げることが考えられるが、例えば−８００ＨＵ程度のＧＧＯ領域が未検出となる場合があり、単純に閾値の調整のみで対応することは困難である。なお、図３は、従来技術に係る課題を説明するための図である。
【００２３】
そこで、本願に係る画像処理装置１００は、以下、詳細に説明する構成により、ボリュームデータからの領域抽出の精度を向上させることを可能にする。具体的には、画像処理装置１００は、前景領域の抽出精度を向上させることにより、ボリュームデータからの領域抽出の精度を向上させる。図４は、第１の実施形態に係る画像処理装置１００の一例を示す図である。図４に示すように、画像処理装置１００は、入力部１１０と、表示部１２０と、記憶部１３０と、制御部１４０とを有する。例えば、画像処理装置１００は、ワークステーションや、任意のパーソナルコンピュータなどであり、医用画像診断装置２００や、画像保管装置３００などとネットワークを介して接続される。
【００２４】
入力部１１０は、マウス等のポインティングデバイス、キーボード等の入力デバイスであり、画像処理装置１００に対する各種操作の入力を操作者から受け付ける。例えば、入力部１１０は、ボリュームデータに含まれる抽出対象の指定操作などを受付ける。表示部１２０は、液晶ディスプレイ等の表示デバイスであり、各種情報を表示する。具体的には、表示部１２０は、操作者から各種操作を受け付けるためのＧＵＩ（Graphical  User  Interface）や、後述する制御部１４０による処理結果を表示する。
【００２５】
記憶部１３０は、例えば、ＲＡＭ（Random  Access  Memory）、フラッシュメモリ（Flash  Memory）等の半導体メモリ素子、又は、ハードディスク、光ディスク等の記憶装置などであり、後述する制御部１４０によって取得されたボリュームデータなどを記憶する。また、記憶部１３０は、後述する制御部１４０によって用いられる種々の情報を記憶する。また、記憶部１３０は、後述する制御部１４０による処理結果を記憶する。例えば、記憶部１３０は、図４に示すように、画像データ記憶部１３１を有する。画像データ記憶部１３１は、後述する制御部１４０によって取得されたボリュームデータを記憶する。
【００２６】
制御部１４０は、例えば、ＣＰＵ（Central  Processing  Unit）やＭＰＵ（Micro  Processing  Unit）等の電子回路、ＡＳＩＣ（Application  Specific  Integrated  Circuit）やＦＰＧＡ（Field  Programmable  Gate  Array）等の集積回路であり、画像処理装置１００の全体制御を行なう。制御部１４０は、図４に示すように、例えば、画像データ取得部１４１と、第１の抽出部１４２と、第２の抽出部１４３とを有する。
【００２７】
画像データ取得部１４１は、医用画像診断装置２００によって生成されたボリュームデータを取得する。例えば、画像データ取得部１４１は、入力部１１０を介して操作者によって指定されたボリュームデータを医用画像診断装置２００又は画像保管装置３００から取得して、画像データ記憶部１３１に格納する。
【００２８】
第１の抽出部１４２は、ボリュームデータから抽出対象の構造を含む第１の領域を抽出する。具体的には、第１の抽出部１４２は、ボリュームデータにおいてボクセルの輝度値が所定の閾値以上である領域を第１の領域として抽出する。
【００２９】
第２の抽出部１４３は、第１の領域をボクセルの特徴量の類似度及び位置関係に基づいて複数の小領域に分割し、ボリュームデータにおける第１の領域以外の領域である第２の領域との類似度及び位置関係に基づいて、複数の小領域から構造に含まれる小領域を抽出することで、構造の領域を抽出する。具体的には、第２の抽出部１４３は、第２の領域と所定の位置関係にある小領域のうち、第２の領域との特徴量分布の類似度が所定の閾値以下、又は、特徴量分布の距離が所定の閾値以上である小領域を構造に含まれる小領域として抽出する。ここで、第２の抽出部１４３は、第１の領域に含まれるボクセルを輝度値及び３次元座標を用いてクラスタリングすることにより、第１の領域を複数の小領域に分割する。
【００３０】
図５は、第１の実施形態に係る画像処理装置１００による処理の手順を示すフローチャートである。図５に示すように、第１の実施形態に係る画像処理装置１００は、ボリュームデータから抽出対象を含む前景領域の抽出処理を、ほぼ自動化された流れで行うことができる。具体的には、第１の実施形態に係る画像処理装置１００は、ボリュームデータから第１の領域を抽出して、抽出した第１の領域を複数の小領域に分割する。そして、画像処理装置１００は、ボリュームデータにおいて第１の領域以外の第２の領域との類似度及び位置関係をもとに各小領域が抽出対象の一部であるか否かを判定し、抽出対象の一部であると判定した小領域を前景領域として抽出する。以下、図６Ａ〜図９を併せて参照しながら、第１の実施形態における処理手順を説明する。なお、以下では、ＣＴボリュームデータから結節の領域を抽出する場合を一例に挙げて説明する。
【００３１】
ステップＳ１０１：まず、画像データ取得部１４１が、操作者による３次元の画像データ（ボリュームデータ）の指定を入力部１１０を介してＧＵＩ上で受け付け、受け付けた指定に従ってボリュームデータを取得する。
【００３２】
ステップＳ１０２：続いて、第１の抽出部１４２が、第１の領域を抽出するための閾値を決定する。図６Ａは、第１の実施形態に係る第１の抽出部１４２による第１の領域の抽出の一例を示す図である。例えば、第１の抽出部１４２は、図６Ａに示すように、画像データ取得部１４１の処理によって取得されたＣＴボリュームデータから結節の領域を含む第１の領域を輝度値によって抽出するための閾値を決定する。ここで、第１の抽出部１４２によって抽出された第１の領域以外の領域を図６Ａに示すように第２の領域とする。なお、第１の抽出部１４２によって実行される輝度値による第１の領域の抽出は、従来技術における前景領域の抽出に相当する。
【００３３】
ここで、第１の抽出部１４２による閾値の決定は、入力部１１０を介して操作者から指定された閾値を、第１の領域を抽出するための閾値として決定する場合であってもよいが、入力された画像ごとに適応的に決定する場合であってもよい。なお、操作者によって指定される閾値としては、例えば、ＧＧＯ領域が可能となるように「−８００ＨＵ」程度の固定を用いる場合であってもよい。
【００３４】
以下、適応的に決定する場合について説明する。図６Ｂは、第１の実施形態に係る第１の抽出部１４２による閾値の決定の一例を説明するための図である。例えば、第１の抽出部１４２は、図６Ｂの（Ａ）に示すように、ＣＴボリュームデータの各ボクセルの輝度値のヒストグラムを生成して、生成したヒストグラムに基づいて閾値を決定する。ここで、ヒストグラム生成に用いる領域は、ボリュームデータ全体であってもよく、或いは、特定の領域に限定したＶＯＩであってもよい。
【００３５】
ここで、Ｓｏｌｉｄ領域やＧＧＯ領域などの解剖学的構造は、それぞれ特有の輝度値を有することから、各構造の輝度値の分布はそれぞれに対応した平均値などのパラメータを持つ確率分布で表現できると仮定すると、図６Ｂの（Ａ）に示すＣＴボリュームデータから生成したヒストグラムは、これら複数の構造に対応する確率分布の和として表現することができる。確率分布は、例えばガンマ分布やｓｔｕｄｅｎｔ’ｓ−ｔ分布などを用いても良いが、以下では、ガウス分布を用いた方法を例に挙げて説明する。例えば、各ボクセルのＣＴボリュームデータにおける座標「ｒ（位置ベクトルｒ）＝（ｘ，ｙ，ｚ）」の輝度値「ｖ＝Ｉ（ｒ）」を確率変数とした混合ガウス分布は、下記の式（１）で表すことができる。
【００３６】
【数１】
0006415878.tif 000002
【００３７】
ここで、式（１）における「ｑ」は確率密度関数を示し、「Ｎ（ｖ｜μｋ，σｋ）」は一次元ガウス分布を示し、「μｋ」は平均値を示し、「σｋ」は標準偏差を示し、「πｋ」は混合率を示し、「Ｋ」は混合要素数を示す。パラメータπｋ，μｋ，σｋは、公知の非特許文献である「Ａ．Ｐ．  Ｄｅｍｐｓｔｅｒ，  Ｎ．Ｍ．  Ｌａｉｒｄ，  ａｎｄ  Ｄ．Ｂ．  Ｒｕｂｉｎ，  “Ｍａｘｉｍｕｍ  Ｌｉｋｅｌｉｈｏｏｄ  ｆｒｏｍ  Ｉｎｃｏｍｐｌｅｔｅ  Ｄａｔａ  ｖｉａ  ｔｈｅ  ＥＭ  Ａｌｇｏｒｉｔｈｍ，”  Ｊｏｕｒｎａｌ  ｏｆ  ｔｈｅ  Ｒｏｙａｌ  Ｓｔａｔｉｓｔｉｃａｌ  Ｓｏｃｉｅｔｙ  Ｓｅｒｉｅｓ  Ｂ，  Ｖｏｌ．３９，  ｐｐ．１−３８，  １９７７．」に開示されている方法により推定することができる。
【００３８】
例えば、各要素（Ｓｏｌｉｄ領域や、ＧＧＯ領域など）のガウス分布（正規分布）の和がヒストグラムを近似するようにパラメータπｋ，μｋ，σｋが推定される。ここで、混合要素数Ｋは任意に設定することができるが、例えば、図６Ｂの（Ｂ）に示すように、「Ｋ＝３」として推定する場合には、各要素それぞれのガウス分布の和が図６Ｂの（Ａ）のヒストグラムを近似するように、各要素のパラメータπｋ，μｋ，σｋが推定される。すなわち、図６Ｂの（Ｂ）に示すように、空気を含む要素（ｋ＝１に対応する要素）のガウス分布と、ＧＧＯ領域（すりガラス状陰影）を含む要素（ｋ＝２に対応する要素）のガウス分布と、血管やＳｏｌｉｄ領域（充実性陰影）を含む要素（ｋ＝３に対応する要素）のガウス分布が得られる。
【００３９】
ここで、抽出対象の構造をＧＧＯ領域、Ｓｏｌｉｄ領域、血管とした場合には、ｋ＝２に対応する要素と、ｋ＝３に対応する要素のみを選択すればよい。そこで、各ボクセルの座標ｒの輝度値ｖが抽出対象の構造である確率ｌ（ｖ）をベイズの定理に基づく下記の式（２）で算出することができる。
【００４０】
【数２】
0006415878.tif 000003
【００４１】
なお、式（２）における分母は全要素（ｋ＝１〜３）の確率分布の混合確率分布を示し、分子はｋ＝２に対応する確率分布とｋ＝３に対応する確率分布の混合確率分布を示す。ここで、ＣＴボリュームデータから第１の領域を抽出するための閾値「ｖth」は、「ｌ（ｖth）＝β」となるように設定する。なお、確率ｌ（ｖ）は、０≦ｌ（ｖ）≦１の範囲をとり事後確率を表現することから、閾値を示すβはベイズ理論における最小誤り率の観点からβ＝０．５としても良く、或いは、０≦β≦１の範囲内の他の値を選択しても良い。
【００４２】
ステップＳ１０３：第１の抽出部１４２は、上述したように閾値を決定することで、例えば、図６Ｂの（Ｂ）に示すように、空気を含む要素のガウス分布とＧＧＯ領域を含むガウス分布との交点付近を閾値「ｖth」とし、第１の領域を「ｖ＞ｖth」を満たすボクセルの集合として抽出する。すなわち、第１の抽出部１４２は、第１の領域を下記の式（３）として抽出する。一方、第１の抽出部１４２は、第２の領域を下記の式（４）として抽出する。なお、式（３）における「Ｆ」は第１の領域を示し、式（４）における「Ｂ」は第２の領域を示す。
【００４３】
【数３】
0006415878.tif 000004
【００４４】
【数４】
0006415878.tif 000005
【００４５】
なお、図６Ａに示す右側の画像は、第１の領域に含まれるボクセルと第２の領域に含まれるボクセルとを二値化した値で表現した画像であり、第１領域に含まれるボクセルを「１」の「白」の階調で表示し、第２の領域に含まれるボクセルを「０」の「黒」の階調で表示したものである。すなわち、図６Ａの例では、第１の抽出部１４２は、「白」の階調で表示された領域を、上述したように、入力されたボリュームデータごとに適応的に閾値を決定することで、第１の領域として抽出する。
【００４６】
ステップＳ１０４：図５に戻って、次に、第２の抽出部１４３が、ステップＳ１０３によって抽出された第１の領域を小領域に分割する。具体的には、第２の抽出部１４３は、第１の領域に含まれるボクセルの特徴量（例えば、輝度値）と３次元座標を用いたクラスタリングにより、第１の領域を複数の小領域に分割する。これにより、特徴量（輝度値）の空間的な連続性が考慮されるため、Ｓｏｌｉｄ領域やＧＧＯ領域などの解剖学的構造の境界に沿った領域分割を行うことができる。
【００４７】
図７は、第１の実施形態に係る第２の抽出部１４３による処理の概念を説明するための模式図である。図７においては、横軸にボリュームデータにおける座標空間を示し、縦軸にボクセルの特徴量（輝度値）を示す。例えば、第２の抽出部１４３は、図７の（Ａ）に示す所定の空間に含まれる複数のボクセルについて、特徴量（輝度値）の分布から１つの小領域としてクラスタリングする。一例を挙げると、第２の抽出部１４３は、図７の（Ｂ）に示すように、ある閾値以下の特徴量を有するボクセルＰ１を所定の空間内（距離内）にある他のボクセルの特徴量（輝度値）の分布に合せて同一の小領域としてクラスタリングする。すなわち、第２の抽出部１４３は、ボクセルＰ１を特徴量（輝度値）の閾値を超えた他のボクセルと同一の小領域として分割する。換言すると第２の抽出部１４３は、ボリュームデータ（第１の領域）に含まれる各ボクセルについて、空間的に近い他のボクセルとの類似度に基づいてどの領域に属するかを判定し、判定結果に応じて領域を区分けする。
【００４８】
ここで、第２の抽出部１４３は、上述したクラスタリングに、例えば、Ｍｅａｎｓｈｉｆｔ法、Ｍｅｄｏｉｄｓｈｉｆｔ法、ＳＬＩＣ法、Ｑｕｉｃｋｓｈｉｆｔ法などを用いることができる。以下、公知の非特許文献である「Ａ．  Ｖｅｄａｌｄｉ  ａｎｄ  Ｓ．  Ｓｏａｔｔｏ，  “Ｑｕｉｃｋ  Ｓｈｉｆｔ  ａｎｄ  Ｋｅｒｎｅｌ  Ｍｅｔｈｏｄｓ  ｆｏｒ  Ｍｏｄｅ  Ｓｅｅｋｉｎｇ，“  ｉｎ  Ｐｒｏｃｅｅｄｉｎｇｓ  ｏｆ  ｔｈｅ  Ｅｕｒｏｐｅａｎ  Ｃｏｎｆｅｒｅｎｃｅ  ｏｎ  Ｃｏｍｐｕｔｅｒ  Ｖｉｓｉｏｎ  （ＥＣＣＶ），  ２００８」に開示されているＱｕｉｃｋｓｈｉｆｔ法を用いる場合について説明する。
【００４９】
Ｑｕｉｃｋｓｈｉｆｔ法はカーネル密度関数の極値として定義されるクラスタ中心を探索する手法であり、（１）カーネル密度関数の計算、（２）最近傍探索、（３）クラスタリングの３ステップで実行される。
【００５０】
（カーネル密度関数の計算）
まず、第２の抽出部１４３は、ボリュームデータに含まれる各ボクセルについて、カーネル密度関数の値を算出する。ここで、ボリュームデータに対してＱｕｉｃｋｓｈｉｆｔ法を適用する場合には、１次元の輝度値ｖ＝Ｉ（ｒ）と３次元画像の座標ｒ＝（ｘ，ｙ，ｚ）を結合した４次元の結合空間Ｒ4を特徴空間とし、輝度値の空間（輝度空間）及び画像空間それぞれにおけるカーネルの積で定義された下記の式（５）によりカーネル密度関数が算出される。
【００５１】
【数５】
0006415878.tif 000006
【００５２】
なお、式（５）における「ｐ（ｖ，ｒ）」がカーネル密度関数を示し、「ｄv」及び「ｈv」がそれぞれ輝度空間の距離尺度とカーネル幅を示し、「ｄr」及び「ｈr」がそれぞれ画像空間の距離尺度とカーネル幅を示し、「ｎ」がボクセル数を示し、「ｋ」がカーネル関数を示し、「Ｃ」が定数を示す。ここで、各空間における距離尺度は、下記の式（６）で示すユークリッド距離によって算出され、カーネル関数ｋは、下記の式（７）に示すEpanechnikovカーネルを用いることが望ましい。これにより、計算量を低減することができ、装置の処理負荷を低減することができる。
【００５３】
【数６】
0006415878.tif 000007
【００５４】
【数７】
0006415878.tif 000008
【００５５】
すなわち、第２の抽出部１４３は、式（５）に示すように、ボリュームデータに含まれる各ボクセルについて、輝度値及び座標それぞれで所定の範囲における類似度を定義したカーネル密度関数の値を算出する。ここで、輝度値及び座標それぞれにおける範囲（カーネル幅）は、抽出対象や撮像条件によって任意に設定することができる。例えば、輝度空間のカーネル幅「ｈv」は、分離する対象の輝度値の差分が設定される。一例を挙げると、空気を含む要素とＧＧＯ領域を含む要素とを精度よく分離するためには、輝度空間のカーネル幅「ｈv」は、空気を示すＣＴ値「−１０００ＨＵ」とＧＧＯ領域を示すＣＴ値「−８００」との差分である「２００」前後に設定することが望ましい。
【００５６】
また、例えば、画像空間のカーネル幅「ｈr」は、抽出対象のサイズに基づいて設定される。一例を挙げると、画像収集時の画素ピッチが「０．６ｍｍ」のボリュームデータから、抽出対象の最小のサイズが「半径：１．５ｍｍ」程度の領域を抽出する場合、画像空間のカーネル幅「ｈr」は、「２．５＝（１．５／０．６）」前後に設定することが望ましい。なお、これらのカーネル幅は、抽出対象の情報や、ボリュームデータの付帯情報（例えば、ＤＩＣＯＭのタグ）に基づいて、第２の抽出部１４３が自動で設定することも可能である。かかる場合には、例えば、第２の抽出部１４３は、予め記憶部１３０に記憶された抽出対象ごとの輝度空間のカーネル幅を読み出すと共に、ボリュームデータから取得した画素ピッチの情報と予め記憶部１３０に記憶された抽出対象のサイズとから画像空間のカーネル幅を算出する。
【００５７】
（最近傍探索）
上述したように、ボリュームデータに含まれる各ボクセルのカーネル密度関数の値を算出すると、第２の抽出部１４３は、上述した式（３）で示される第１の領域に含まれる各ボクセル「ｉ」の近傍「Ｄi」で確率密度が最も増加するボクセル「Ｎ（ｉ）」を下記の式（８）により探索する。
【００５８】
【数８】
0006415878.tif 000009
【００５９】
ここで、近傍Ｄiは任意に設定することができるが、探索に係る処理負荷の低減と、離れたボクセルが割り当てられることを抑止するために、下記の式（９）で定義するように、上述したカーネル幅の範囲を用いることが望ましい。
【００６０】
【数９】
0006415878.tif 000010
【００６１】
ここで、式（８）及び（９）に示す「ｊ」は、ボクセル「ｉ」の近傍「Ｄi」に含まれるボクセルを示し、式（９）に示すように、ボクセル「ｉ」との輝度値の差が輝度空間のカーネル幅以内であり、かつ、ボクセル「ｉ」との座標間の距離が画像空間のカーネル幅以内のボクセルである。このとき、ボクセル「ｊ」は、上述した式（４）で示す第２の領域に含まれるボクセルであってもよい。
【００６２】
また、式（８）に示す「ｄ（ｊ||ｉ）」は、ボクセル「ｉ」とボクセル「ｊ」との距離尺度であり、下記の式（１０）に示す輝度空間と画像空間におけるユークリッド距離の重み付和によって算出される。
【００６３】
【数１０】
0006415878.tif 000011
【００６４】
上述したように、第２の抽出部１４３は、第１の領域に含まれる各ボクセル「ｉ」の近傍「Ｄi」で確率密度が最も増加するボクセル「Ｎ（ｉ）」を探索する。すなわち、第２の抽出部１４３は、ボクセル「ｉ」の近傍で類似したボクセルの密度が最も高い位置に相当するボクセル「Ｎ（ｉ）」を探索する。
【００６５】
（クラスタリング）
上述した最近傍探索を行うと、第２の抽出部１４３は、次に最近傍探索の結果を用いて第１の領域を複数の小領域に分割する。図８Ａは、第１の実施形態に係る第２の抽出部１４３による領域分割の一例を示す図である。例えば、第２の抽出部１４３は、図８Ａに示すように、第１の抽出部１４２によって抽出された第１の領域を複数の小領域に分割する。ここで、第２の抽出部１４３は、最近傍探索の結果を用いたボクセルの遷移操作を再帰的に繰り返すことで各小領域の中心（クラスタ中心）を探索するとともに、各小領域に含まれるボクセルを抽出する。
【００６６】
図８Ｂは、第１の実施形態に係る第２の抽出部１４３によるクラスタリングの一例を説明するための模式図である。例えば、第２の抽出部１４３は、図８Ｂに示すように、第１の領域に含まれるボクセル「ｉ」をグラフ上の１つのノードとして、ボクセル「ｉ」と最近傍探索によって探索されたＮ（ｉ）との間にエッジを張る。第２の抽出部１４３は、第１の領域に含まれるすべてのボクセルについて同様にエッジを張ることで画像全体を複数の木構造で表現する。ここで、クラスタ中心は、確率密度関数の極大値として定義されるため、その近傍領域により確率密度が高い点が存在しない点であると解釈できる。すなわち、クラスタ中心の点に相当するボクセルは、自身の点がＮ（ｉ）となるためエッジが張られず、木構造の根に相当することになる。
【００６７】
そこで、第２の抽出部１４３は、第１の領域に含まれるすべてのボクセル「ｉ」に対して、図８Ｂに示す「ｉ→Ｎ（ｉ）→Ｎ（Ｎ（ｉ））」のような遷移操作を再帰的に繰り返し実行することで木構造の根「Ｒ（ｉ）」の探索を行う。すなわち、第２の抽出部１４３は、図８Ｂに示すように、ボクセル「ｉ」をスタートとして、最近傍探索の結果であるボクセル「Ｎ（ｉ）」に遷移する。そして、第２の抽出部１４３は、ボクセル「Ｎ（ｉ）」の最近傍探索の結果であるボクセル「Ｎ（Ｎ（ｉ））」に遷移する。そして、第２の抽出部１４３は、最近傍探索の結果が自身のボクセルになる木構造の根「Ｒ（ｉ）」を探索する。
【００６８】
ここで、根のボクセル位置「ｒR(i)」が第２の領域に含まれる（ｒR(i)∈Ｂ）場合、第２の抽出部１４３は、スタートのボクセル「ｉ」を第１の領域から除外する。一方、根のボクセル位置「ｒR(i)」が第１の領域に含まれる（ｒR(i)∈Ｆ）場合、第２の抽出部１４３は、スタートのボクセル「ｉ」に対して（根）「Ｒ（ｉ）」をクラスタ中心としたクラスタ（小領域）を割り当てる。このとき、ボクセル「ｉ」の輝度値「ｖi」と根のボクセル「Ｒ（ｉ）」の輝度値「ｖR(i)」の差が下記の式（１１）となる場合には、木を分割することで新しいクラスタとする。式（１１）において、「ａ」は予め設定される定数であり、例えば、「ａ＝１」や、「ａ＝２」のように設定される。これにより、解剖学的構造の境界での領域分割をより正確に行うことを可能にする。
【００６９】
【数１１】
0006415878.tif 000012
【００７０】
第２の抽出部１４３は、上述したクラスタリングを第１の領域に含まれるすべてのボクセルに対して実行することで、第１の領域をグラフ上の根「Ｒ（ｉ）」それぞれをクラスタ中心とした複数の小領域に分割する。
【００７１】
ステップＳ１０５：図５に戻って、第１の領域を小領域に分割すると、第２の抽出部１４３は、特定の情報を受け付けたか否かの判定を行う。具体的には、第２の抽出部１４３は、ボリュームデータに含まれる１つ又は複数のボクセルが抽出対象の構造に該当するか否かを特定するための特定情報（位置情報）を受け付けたか否かを判定する。ここで、位置情報は、例えば、操作者のマウスクリックやドラッグ操作などによって指定された点又は領域であってもよく、或いは、ＣＡＤによって自動検出された点又は領域であってもよい。
【００７２】
ここで、特定情報を受け付けた場合には（ステップＳ１０５、Ｙｅｓ）、第２の抽出部１４３は、特定情報を参照しつつ、対象構造の一部である小領域を抽出する（ステップＳ１０６）。すなわち、第２の抽出部１４３は、特定情報を受け付けたボクセル又は領域を含む小領域を前景領域（対象構造の一部）として抽出する。一方、特定情報を受け付けていない場合には（ステップＳ１０５、Ｎｏ）、対象構造の一部である小領域を抽出する（ステップＳ１０７）。
【００７３】
ステップＳ１０６、ステップＳ１０７：ステップＳ１０５において特定情報を受け付けていない場合、及び、処理対象の小領域が特定情報を受け付けたボクセル又は領域を含んでいない場合には、第２の抽出部１４３は、分割した小領域が抽出対象の構造の一部であるか否かを判定して、抽出対象の構造の一部と判定した小領域を前景領域として抽出する。図９は、第１の実施形態に係る第２の抽出部１４３による抽出処理の一例を示す図である。例えば、第２の抽出部１４３は、図９に示すように、分割した小領域のうち、抽出対象の構造の一部ではない小領域を除外した小領域を前景領域として抽出する。
【００７４】
ここで、小領域に対する除外処理は、第１の領域を小領域に分割する際に用いたＱｕｉｃｋｓｈｉｆｔ法に基づく方法により実行される。なお、小領域に対する除外処理は、公知の非特許文献である「Ｊ．  Ｎｉｎｇ，  Ｌ．  Ｚｈａｎｇ，  Ｄ．  Ｚｈａｎｇ，  ａｎｄ  Ｃ．  Ｗｕ，  “Ｉｎｔｅｒａｃｔｉｖｅ  ｉｍａｇｅ  ｓｅｇｍｅｎｔａｔｉｏｎ  ｂｙ  ｍａｘｉｍａｌ  ｓｉｍｉｌａｒｉｔｙ  ｂａｓｅｄ  ｒｅｇｉｏｎ  ｍｅｒｇｉｎｇ，”  Ｐａｔｔｅｒｎ  Ｒｅｃｏｇｎｉｔｉｏｎ，  ｖｏｌ．  ４３，  Ｎｏ．  ２，  ｐｐ．  ４４５−４５６，  ２０１０．」に開示されている方法でもＱｕｉｃｋｓｈｉｆｔ法と同様に実行することができる。
【００７５】
Ｑｕｉｃｋｓｈｉｆｔ法はグラフ構造の各ノードにおいてカーネル密度関数の値を算出し、密度の増加する方向にグラフ上を遷移しながら極大値の探索を行う手法であると解釈できる。上述した小領域の分割においては、各ノードがボクセルに対応し、ノード間の距離を輝度値及び画像上の距離として定義した。小領域の除外処理においては、各ノードを小領域としたうえでその位置関係（隣接関係）からグラフ構造を構成し、ノード（領域）間の距離を適切に定義すればよい。基本的な流れは、上述した小領域の分割と同様に、（１）カーネル密度関数の計算、（２）最近傍探索、（３）クラスタリングの３ステップで実行される。
【００７６】
（カーネル密度関数の計算）
まず、第２の抽出部１４３は、得られたＭ個の小領域に第２の領域を加えたＭ＋１個の領域sj（j=１,２,・・・，Ｍ＋１)ごとに、下記の式（１２）で定義されるカーネル密度関数の値を算出する。
【００７７】
【数１２】
0006415878.tif 000013
【００７８】
ここで、式（１２）における「ｈ（ｓ）」は小領域「ｓ」に含まれるボクセルの輝度値のヒストグラムを示し、「ｄH」はヒストグラム間の距離を示し，「ｄG」は小領域間のグラフ距離を定義する関数を示し、「ｈH」及び「ｈG」はそれぞれヒストグラム間のカーネル幅及び小領域間のカーネル幅を示す。ｄHにはKullback-Leibler距離などを用いても良いが、本実施形態では下記の式（１３）で定義されるHellinger距離を用いる場合について説明する。
【００７９】
【数１３】
0006415878.tif 000014
【００８０】
ここで、式（１３）における「Ｎb」はヒストグラムのビン数を示し、「Ｂ」はBattacharyya係数を示す。ｄGには小領域内に含まれるボクセルの重心座標間のユークリッド距離を用いても良いが、領域が凸でない場合（領域にへこみ部分がある場合や、中心付近に穴（他の小領域）がある場合など）に重心の座標が領域外となってしまう場合がある。例えば、領域の中心付近に他の小領域がある場合、領域の重心が他の小領域内になってしまう場合がある。
【００８１】
そこで、小領域の隣接関係を表す隣接行列（各ｉ、ｊ成分は領域ｓiと領域ｓjが隣接している場合に１、隣接していない場合に０となる）に基づきグラフ構造を生成し，ｄGをグラフ上の距離、例えば最小全域木（Minimum  Spanning  Tree）や最短経路（２つのノードをつなぐ経路上で通過する最小のエッジ数）を基に設定すると良い。すなわち、隣接するノード（小領域）間でのみエッジを張ったグラフ構造に基づいて、グラフ上の距離を設定する。本実施形態では、最短経路を用いる場合について説明する。ここで、カーネル関数をEpanechnikovカーネルとし、処理負荷を低減するためにカーネル幅をｈH=１、ｈG=２とすると、カーネル密度関数は下記の式（１４）と定義することができる。
【００８２】
【数１４】
0006415878.tif 000015
【００８３】
ここで、Ｓ：＝｛ｓi|ｄG（ｓ，ｓi)＜２}である。また、ｄG（ｓ，ｓi）はエッジ数であり整数値に限定されることから、ｄG（ｓ，ｓi）＜２となるのはｄG（ｓ，ｓi）＝０と、ｄG（ｓ，ｓi）＝１であることに注意すると、カーネル密度関数は、さらに下記の式（１５）と定義することができる。なお、式（１５）においては、共通の定数部分は省略されている。
【００８４】
【数１５】
0006415878.tif 000016
【００８５】
ここで、式（１５）における「Ｎs」は小領域ｓに隣接する小領域の集合を示す。このように設定することで、カーネル密度関数は、式（１５）に示すように、各領域ｓに隣接する領域とのBattacharyya係数の総和として算出することができる。
【００８６】
（最近傍探索）
上述したように、各小領域のカーネル密度関数の値を算出すると、第２の抽出部１４３は、各小領域「ｓj」の近傍「Ｄsj」で確率密度が最も増加する小領域「Ｎ（ｓj）」を下記の式（１６）により探索する。
【００８７】
【数１６】
0006415878.tif 000017
【００８８】
ここで、近傍Ｄsjは任意に設定することができるが、探索に係る処理負荷の低減と、離れた小領域が割り当てられることを抑止するために、下記の式（１７）で定義するように、設定することが望ましい。
【００８９】
【数１７】
0006415878.tif 000018
【００９０】
ここで、式（１７）に示す「ＴH」はヒストグラム間の距離に関する閾値を示し、「ＴG」はグラフ距離に関する閾値を示す。これらの閾値は、それぞれHellinger距離と最短経路とを用いる場合に、下記の式（１８）として設定することも可能である。
【００９１】
【数１８】
0006415878.tif 000019
【００９２】
上述したように、第２の抽出部１４３は、各小領域「ｓj」の近傍「Ｄsj」で確率密度が最も増加する小領域「Ｎ（sj）」を探索する。すなわち、第２の抽出部１４３は、小領域「ｓj」の近傍で類似した小領域の密度が最も高い位置に相当する小領域「Ｎ（sj）」を探索する。
【００９３】
（クラスタリング）
上述した最近傍探索を行うと、第２の抽出部１４３は、次に最近傍探索の結果を用いて前景領域を抽出する。ここで、第２の抽出部１４３は、最近傍探索の結果を用いた小領域の遷移操作を再帰的に繰り返すことでクラスタ中心を探索するとともに、各小領域が抽出対象の構造の一部であるか否か判定して、一部である場合に前景領域として抽出する（一部ではない場合に、第１の領域から除外する）。
【００９４】
例えば、第２の抽出部１４３は、各小領域「ｓj」をグラフ上の１つのノードとして、小領域「ｓj」と最近傍探索によって探索されたＮ（ｓj）との間にエッジを張る。第２の抽出部１４３は、すべての小領域について同様にエッジを張ることで画像全体を複数の木構造で表現する。ここで、クラスタ中心は、確率密度関数の極大値として定義されるため、その近傍により確率密度が高い小領域が存在しない小領域であると解釈できる。すなわち、クラスタ中心の小領域は、自身がＮ（ｓj）となるためエッジが張られず、木構造の根に相当することになる。
【００９５】
そこで、第２の抽出部１４３は、すべての小領域「ｓj」に対して、「ｓj→Ｎ（ｓj）→Ｎ（Ｎ（ｓj））」のような遷移操作を再帰的に繰り返し実行することで木構造の根「Ｒ（ｓj）」の探索を行う。すなわち、第２の抽出部１４３は、小領域「ｓj」をスタートとして、最近傍探索の結果である小領域「Ｎ（ｓj）」に遷移する。そして、第２の抽出部１４３は、小領域「Ｎ（ｓj）」の最近傍探索の結果であるボクセル「Ｎ（Ｎ（ｓj））」に遷移する。そして、第２の抽出部１４３は、最近傍探索の結果が自身の領域になる木構造の根「Ｒ（ｓj）」を探索する。
【００９６】
ここで、根の小領域「ｓR(ｓj)」が第２の領域の場合、第２の抽出部１４３は、スタートの小領域「ｓj」を第１の領域から除外する。一方、根の小領域「ｓR(ｓj)」が第１の領域の場合、第２の抽出部１４３は、スタートの小領域「ｓj」を前景領域として抽出する。第２の抽出部１４３は、上述したクラスタリングをすべての小領域に対して実行することで、前景領域を抽出する。ここで、小領域の内の一つ又は複数が、抽出対象の構造に該当するための特定情報（位置情報）に対応していた場合には、特定情報に応じた抽出処理が実行される。例えば、操作者により手動、又はＣＡＤによって自動的に指定された抽出対象構造のボクセル位置を含む小領域は、たとえ根の領域が第２領域であったとしても第１領域から除外しない、というようにすることもできる。
【００９７】
上述したように、第１の実施形態によれば、第１の抽出部１４２は、３次元画像データから抽出対象の構造を含む第１の領域を抽出する。第２の抽出部１４３は、第１の領域をボクセルの特徴量の類似度及び位置関係に基づいて複数の小領域に分割し、３次元画像データにおける第１の領域以外の領域である第２の領域との類似度及び位置関係に基づいて、複数の小領域から構造に含まれる小領域を抽出することで、構造の領域を抽出する。従って、第１の実施形態に係る画像処理装置１００は、ボクセル毎の輝度値だけでなく、領域同士の輝度値の分布の類似度と空間的隣接関係を考慮したことで解剖学的構造の境界に沿った領域抽出ができ、３次元画像データからの領域抽出の精度を向上させることを可能にする。
【００９８】
（第２の実施形態）
第１の実施形態では、抽出対象の構造の一部であると判定した小領域を前景領域として抽出する場合について説明した。第２の実施形態では、抽出した小領域がいずれの抽出対象構造であるかを識別する場合について説明する。図１０は、第２の実施形態に係る画像処理装置１００の一例を示す図である。図１０に示すように、第２の実施形態に係る画像処理装置１００は、記憶部１３０に推定基準記憶部１３２を新たに有する点と、制御部１４０に識別部１４４を新たに有する点とが第１の実施形態と異なる。以下、これらを中心に説明する。
【００９９】
推定基準記憶部１３２は、抽出した小領域がいずれの抽出対象構造であるかを識別するための推定基準を記憶する。具体的には、推定基準記憶部１３２は、種々の手法によって生成された分類器を記憶する。例えば、抽出対象がＳｏｌｉｄ領域、ＧＧＯ領域、血管領域である場合、輝度値や領域内のテクスチャ、形状が異なる為、特徴量としてそれらを判別するための分類器が記憶される。
【０１００】
Ｓｏｌｉｄ領域はＧＧＯ領域と比較して高輝度かつ内部のテクスチャが平坦であることから、領域内の輝度値の分布、勾配の情報を考慮した特徴量を用いることで識別可能である。血管領域はＳｏｌｉｄ領域と同程度の輝度値を有するが、管状の形状を有し、塊形状の結節領域とは異なるため、そのような３次元形状を捉える特徴量を用いることで識別可能である。一例を挙げると、推定基準記憶部１３２は、領域内の輝度値の分布、ヘッセ行列の固有値の平均値や分散などの統計量、Ｔｏｔａｌ  ｖａｒｉａｔｉｏｎ、Ｌｏｃａｌ  Ｂｉｎａｒｙ  Ｐａｔｔｅｒｎ、Ｇｒａｙ−Ｌｅｖｅｌ  Ｃｏ−ｏｃｃｕｒｒｅｎｃｅ  Ｍａｔｒｉｘ、Ｈｉｓｔｏｇｒａｍ  ｏｆ  Ｏｒｉｅｎｔｅｄ  Ｇｒａｄｉｅｎｔ、Ｈａａｒ−ｌｉｋｅ特徴、Ｈｉｇｈｅｒ−ｏｒｄｅｒ  Ｌｏｃａｌ  Ａｕｔｏ  Ｃｏｒｒｅｌａｔｉｏｎなどを特徴量として判別する分類器を記憶する。また、推定基準記憶部１３２は、上述した特徴量の組み合わせや、Ｂａｇ  ｏｆ  ｆｅａｔｕｒｅやＳｐａｒｓｅ  ｃｏｄｉｎｇを用いて学習したテクスチャパターンの辞書、又はそれを利用した特徴量を判別する分類器を記憶する。
【０１０１】
かかる場合には、例えば、制御部１４０は、まず、描出された構造がいずれの構造であるかが判別された複数の教師有画像における各構造の特徴量をそれぞれ抽出し、抽出した特徴量を用いて分類器を生成して推定基準記憶部１３２に格納する。分類器は識別部１４４によって識別されるごとに更新されても良い。
【０１０２】
識別部１４４は、小領域ごとに上述した種々の特徴量を抽出し、推定基準記憶部１３２によって記憶された分類器を用いて小領域がいずれの抽出対象構造であるかを識別する。また、識別部１４４は、小領域を識別した後に、いずれの抽出対象構造であるかの判別結果を、入力部１１０を介して医師などから受け付けても良い。そして、識別部１４４は、分類器による識別結果と、受け付けた判別結果とを用いて、推定基準記憶部１３２に格納した分類器を更新することもできる。このように、識別部１４４が小領域の特徴量を識別するごとに分類器を更新することで、より精度の高い分類器を生成することができる。
【０１０３】
なお、分類器の生成には種々の手法を用いることができ、例えばＤｉｓｃｒｉｍｉｎａｎｔ  Ａｎａｌｙｓｉｓ、Ｌｏｇｉｓｔｉｃ  Ｒｅｇｒｅｓｓｉｏｎ、Ｓｕｐｐｏｒｔ  Ｖｅｃｔｏｒ  Ｍａｃｈｉｎｅ、Ｎｅｕｒａｌ  Ｎｅｔｗｏｒｋ、Ｒａｎｄｏｍｉｚｅｄ  Ｔｒｅｅ、Ｓｕｂ−Ｓｐａｃｅ  Ｍｅｔｈｏｄなどが挙げられる。ここで、Ｃｏｎｄｉｔｉｏｎａｌ  Ｒａｎｄｏｍ  ＦｉｅｌｄｓやＧｒａｐｈ  Ｃｕｔと組み合わせることで各小領域の隣接関係を考慮した識別を行うと良い。
【０１０４】
図１１は、第２の実施形態に係る画像処理装置による処理の手順を示すフローチャートである。
【０１０５】
ステップＳ２０１：ステップＳ１０１〜ステップＳ１０７にて小領域が抽出されると、識別部１４４は、抽出された小領域がいずれの抽出対象構造であるかを識別する。
【０１０６】
ステップＳ２０２：その後、識別部１４４は、識別結果に基づいて、抽出対象構造ごとに区分けした領域を出力する。図１２は、第２の実施形態に係る識別部１４４による処理の一例を説明するための図である。例えば、識別部１４４は、図１２に示すように、第２の抽出部１４３によって抽出された各小領域について、Ｓｏｌｉｄ領域、ＧＧＯ領域、血管領域のいずれであるかを識別し、識別結果が同一の構造で隣接する小領域を１つの領域として区分けした結果を出力する。すなわち、識別部１４４は、隣接する小領域でＧＧＯ領域と識別した小領域を１つのＧＧＯ領域として出力する。同様に、識別部１４４は、Ｓｏｌｉｄ領域、血管領域についてもそれぞれ区分けした領域を出力する。
【０１０７】
上述したように、第２の実施形態によれば、識別部１４４が、構造の領域に含まれる小領域がそれぞれ抽出対象のいずれであるかを識別する。従って、第２の実施形態に係る画像処理装置１００は、ボクセルごとの識別でなく領域ごとの識別を行うことで、領域内の統計的な特徴量を利用できノイズに対して頑健な識別を行うことを可能にする。また、第２の抽出部１４３で生成された小領域は解剖学的構造の境界に沿って区分されており、図３に示したような過検出領域は除外されている為、特徴抽出の際に不要な情報が混入せず、図１２に示したように精度良く識別を行うことを可能にする。
【０１０８】
（第３の実施形態）
なお、実施形態は、上述した第１及び第２の実施形態に限られるものではない。
【０１０９】
上述した実施形態では、Ｘ線ＣＴ装置によって収集されたＣＴボリュームデータから領域を抽出する場合について説明した。しかしながら、実施形態はこれに限定されるものではなく、ＭＲＩ装置や超音波診断装置などのその他のモダリティによって収集されたボリュームデータから領域を抽出する場合であってもよい。
【０１１０】
また、上述した実施形態では、画像データ取得部１４１が、画像保管装置３００又は医用画像診断装置２００からボリュームデータを取得する場合について説明した。しかしながら、実施形態はこれに限定されるものではなく、例えば、医師がフラッシュメモリや外付けハードディスクなどの可搬性の記憶媒体で医用画像データを持ち運び、画像処理装置１００の画像データ記憶部１３１に格納する場合であってもよい。かかる場合には、画像データ取得部１４１によるボリュームデータの取得は実行されなくてもよい。
【０１１１】
また、上述した実施形態では、画像処理装置１００について説明したが、実施形態はこれに限定されるものではなく、例えば、図４に示す画像処理装置１００の制御部１４０が図１に示す医用画像診断装置２００に組み込まれ、医用画像診断装置にて上述した処理が実行される場合であってもよい。
【０１１２】
また、上述した実施形態の中で示した処理手順に示された指示は、ソフトウェアであるプログラムに基づいて実行されることが可能である。汎用コンピュータが、このプログラムを予め記憶しておき、このプログラムを読み込むことにより、上述した実施形態の画像処理装置１００による効果と同様の効果を得ることも可能である。上述した実施形態で記述された指示は、コンピュータに実行させることのできるプログラムとして、磁気ディスク（フレキシブルディスク、ハードディスクなど）、光ディスク（ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ、ＤＶＤ−ＲＯＭ、ＤＶＤ±Ｒ、ＤＶＤ±ＲＷなど）、半導体メモリ、又はこれに類する記録媒体に記録される。コンピュータ又は組み込みシステムが読み取り可能な記憶媒体であれば、その記憶形式は何れの形態であってもよい。コンピュータは、この記録媒体からプログラムを読み込み、このプログラムに基づいてプログラムに記述されている指示をＣＰＵで実行させれば、上述した実施形態の画像処理装置１００と同様の動作を実現することができる。また、コンピュータがプログラムを取得する場合又は読み込む場合は、ネットワークを通じて取得又は読み込んでもよい。
【０１１３】
また、記憶媒体からコンピュータや組み込みシステムにインストールされたプログラムの指示に基づきコンピュータ上で稼働しているＯＳ（Operating  System）や、データベース管理ソフト、ネットワーク等のＭＷ（Middleware）等が、上述した実施形態を実現するための各処理の一部を実行してもよい。更に、記憶媒体は、コンピュータあるいは組み込みシステムと独立した媒体に限らず、ＬＡＮ（Local  Area  Network）やインターネット等により伝達されたプログラムをダウンロードして記憶又は一時記憶した記憶媒体も含まれる。また、記憶媒体は１つに限られず、複数の媒体から、上述した実施形態における処理が実行される場合も、実施形態における記憶媒体に含まれ、媒体は何れの構成であってもよい。
【０１１４】
なお、実施形態におけるコンピュータ又は組み込みシステムは、記憶媒体に記憶されたプログラムに基づき、上述した実施形態における各処理を実行するためのものであって、パソコン、マイコン等の１つからなる装置、複数の装置がネットワーク接続されたシステム等の何れの構成であってもよい。また、実施形態におけるコンピュータとは、パソコンに限らず、情報処理機器に含まれる演算処理装置、マイコン等も含み、プログラムによって実施形態における機能を実現することが可能な機器、装置を総称している。
【０１１５】
（ハードウェア構成）
図１３は、実施形態に係る画像処理装置のハードウェア構成を示す図である。上述した実施形態に係る画像処理装置は、ＣＰＵ（Central  Processing  Unit）４０等の制御装置と、ＲＯＭ（Read  Only  Memory）５０やＲＡＭ（Random  Access  Memory）６０等の記憶装置と、ネットワークに接続して通信を行う通信Ｉ／Ｆ７０と、各部を接続するバス８０とを備えている。
【０１１６】
上述した実施形態に係る画像処理装置で実行されるプログラムは、ＲＯＭ５０等に予め組み込まれて提供される。また、上述した実施形態に係る画像処理装置で実行されるプログラムは、コンピュータを上述した画像処理装置の各部（例えば、第２の抽出部１４３など）として機能させ得る。このコンピュータは、ＣＰＵ４０がコンピュータ読取可能な記憶媒体からプログラムを主記憶装置上に読み出して実行することができる。
【０１１７】
以上、説明したとおり、第１〜第３の実施形態よれば、３次元画像データからの領域抽出の精度を向上させることができる。
【０１１８】
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【０１１９】
１００  画像処理装置
１４２  第１の抽出部
１４３  第２の抽出部
【図１】
0006415878.tif 000020
【図２】
0006415878.tif 000021
【図３】
0006415878.tif 000022
【図４】
0006415878.tif 000023
【図５】
0006415878.tif 000024
【図６Ａ】
0006415878.tif 000025
【図６Ｂ】
0006415878.tif 000026
【図７】
0006415878.tif 000027
【図８Ａ】
0006415878.tif 000028
【図８Ｂ】
0006415878.tif 000029
【図９】
0006415878.tif 000030
【図１０】
0006415878.tif 000031
【図１１】
0006415878.tif 000032
【図１２】
0006415878.tif 000033
【図１３】
0006415878.tif 000034