(57)【要約】領域分割方法は、空間係数、強度係数及び体積形状係数をクラスタリングするステップを備え、医学的な病巣を自動的に領域分割する。アルゴリズムは以下のステップ、即ち、（１）画像における各ボクセルについて体積形状係数（ＳＩ）を計算するステップ、（２）５次元特徴ベクトルのセットを形成するために、ＳＩ特徴を強度範囲及び空間的位置（ｘ，ｙ，ｚ）と組み合わせるステップ、（３）５次元特徴ベクトルのセットをクラスタに分類するステップ、（４）隣接する領域又はモードを結合するクラスタリングアルゴリズムが提供する強度モードマップに空間特徴だけでなく形状特徴を考慮する改良期待値最大化アルゴリズム（ＥＭ）を用いるステップを有する。結合した空間−強度−形状特徴は、病巣又は腫瘍等の重要な解剖学的構造の領域分割に優良な情報を提供する。
2011526508.tif 000002
【特許請求の範囲】
【請求項１】
３次元医療用画像の領域分割をコンピュータによって実現する方法において、
ａ．前記医療用画像の空間的な変動を示す医療用画像データを受け取るステップと、
ｂ．該医療用画像データから形状データを導出するステップと、
ｃ．複数のモードを生成するために、前記医療用画像データ及び形状データを用いて画像データをクラスタリングするステップと、
ｄ．前記複数のモードから導出され、前記医療用画像に関する領域分割データを出力するステップと
を備えることを特徴とする方法。
【請求項２】
前記医療用画像データは前記医療用画像の空間的な強度の変動を示すこと
を特徴とする請求項１に記載の方法。
【請求項３】
前記形状データは前記医療用画像の体積形状係数を示すこと
を特徴とする請求項１又は請求項２に記載の方法。
【請求項４】
前記形状データはガウス曲率に基づくこと
を特徴とする請求項１から請求項３のいずれか１つに記載の方法。
【請求項５】
前記形状データは平均曲率に基づくこと
を特徴とする請求項１から請求項４のいずれか１つに記載の方法。
【請求項６】
前記形状データはヘッシアン行列から決定されること
を特徴とする請求項１から請求項５のいずれか１つに記載の方法。
【請求項７】
前記形状データはヘッシアン行列の固有値から計算されること
を特徴とする請求項６に記載の方法。
【請求項８】
前記クラスタリングするステップは、夫々が前記医療画像内での空間的位置を構成する一組のベクトルと、前記医療用画像に対応する医療用画像データと、前記医療用画像に対応する形状データをクラスタリングすること
を特徴とする請求項１から請求項７のいずれか１つに記載の方法。
【請求項９】
前記クラスタリングするステップは前記医療用画像データに対応するモードマップを生成すること
を特徴とする請求項１から請求項８のいずれか１つに記載の方法。
【請求項１０】
前記クラスタリングするステップは前記形状データに対応するモードマップを生成すること
を特徴とする請求項１から請求項９のいずれか１つに記載の方法。
【請求項１１】
前記クラスタリングするステップは、平均シフトクラスタリングを行うステップを備えること
を特徴とする請求項１から請求項１０のいずれか１つに記載の方法。
【請求項１２】
前記クラスタリングするステップは、減法クラスタリングを行うステップを備えること
を特徴とする請求項１から請求項１０のいずれか１つに記載の方法。
【請求項１３】
前記クラスタリングするステップは、適合可能なカーネル窓の大きさを用いること
を特徴とする請求項１から請求項１２のいずれか１つに記載の方法。
【請求項１４】
前記モードのうち、隣接するモードを結合するステップを更に備え、
前記領域分割データは前記結合したモードに基づくこと
を特徴とする請求項１から請求項１３のいずれか１つに記載の方法。
【請求項１５】
前記結合するステップは期待値最大化法を行うステップを備えること
を特徴とする請求項１４に記載の方法。
【請求項１６】
前記期待値最大化法を行うステップは前記形状データに基づくこと
を特徴とする請求項１５に記載の方法。
【請求項１７】
前記期待値最大化法を行うステップは空間的な事前確率を計算するステップを有すること
を特徴とする請求項１５又は請求項１６に記載の方法。
【請求項１８】
前記空間的な事前確率は空間的な制限を取り込むこと
を特徴とする請求項１７に記載の方法。
【請求項１９】
前記空間的な制限はマルコフ確率場によって課されること
を特徴とする請求項１８に記載の方法。
【請求項２０】
前記空間的な制限はギブス確率場によって課されること
を特徴とする請求項１８又は請求項１９に記載の方法。
【請求項２１】
前記空間的な事前確率は形状の制限を取り込むこと
を特徴とする請求項１７から請求項２０のいずれか１つに記載の方法。
【請求項２２】
前記モードを結合するステップはグラフ理論ステップを備えること
を特徴とする請求項１４に記載の方法。
【請求項２３】
前記ステップのうち、少なくとも１つは自動的に実施されること
を特徴とする請求項１から請求項２２のいずれか１つに記載の方法。
【請求項２４】
前記医療用画像は、少なくとも結腸の一部の画像であること
を特徴とする請求項１から請求項２３のいずれか１つに記載の方法。
【請求項２５】
前記領域分割データにより前記医療用画像からポリープを抽出すること
を特徴とする請求項２４に記載の方法。
【請求項２６】
前記医療用画像は少なくとも肺の一部の画像であること
を特徴とする請求項１から請求項２３のいずれか１つに記載の方法。
【請求項２７】
前記領域分割データにより前記医療用画像から病巣を抽出すること
を特徴とする請求項２６に記載の方法。
【請求項２８】
前記医療用画像はコンピュータ断層撮影（ＣＴ；Computed Tomography）画像であること
を特徴とする請求項１から請求項２７のいずれか１つに記載の方法。
【請求項２９】
前記医療用画像は磁気共鳴（ＭＲ；Magnetic Resonance）画像であること
を特徴とする請求項１から請求項２７のいずれか１つに記載の方法。
【請求項３０】
前記医療用画像はＸ線画像であること
を特徴とする請求項１から請求項２７のいずれか１つに記載の方法。
【請求項３１】
前記医療用画像は超音波画像であること
を特徴とする請求項１から請求項２７のいずれか１つに記載の方法。
【請求項３２】
３次元医療用画像を領域分割する方法において、
ａ．前記医療用画像の特性の空間的な変動を示す医療用画像データを受け取るステップと、
ｂ．複数のモードを生成するために、前記医療用画像データをクラスタリングするステップと、
ｃ．前記複数のモードのうち、隣接するモードを結合するステップと、
ｄ．前記結合されたモードから導出され、前記医療用画像に関する領域分割データを出力するステップと
を備えることを特徴とする方法。
【請求項３３】
３次元医療用画像を領域分割する方法において、
ａ．前記医療用画像の複数のボクセルの強度を示す医療用画像データを受け取るステップと、
ｂ．該医療用画像データから前記複数のボクセル夫々について形状データを導出するステップと、
ｃ．前記複数のボクセルに対応し、夫々が対応するボクセルについて３次元座標、強度データ及び形状データを有する複数のベクトルを構成するステップと、
ｄ．複数のモードを生成するために前記複数のベクトルをクラスタリングするステップと、
ｅ．前記複数のモードのうち、隣接するモードを結合するステップと、
ｆ．前記結合した複数のモードから導出され、前記医療用画像に関する領域分割データを出力するステップと
を備えることを特徴とする方法。
【請求項３４】
３次元医療用画像を領域分割する方法において、
形状データを含む、前記医療用画像の空間的特徴データをクラスタリングするステップと、
前記クラスタに基づいて前記医療用画像を領域分割するステップと
を備えることを特徴とする方法。
【請求項３５】
適切に準備されたプロセッサによって実行された場合に、請求項１から請求項３４のいずれか１つに記載の方法を実行するために準備されたプログラムコード手段を備えること
を特徴とするコンピュータプログラム。
【請求項３６】
請求項３５に記載のコンピュータプログラムを備えるコンピュータプログラム製品。
【請求項３７】
請求項１から請求項３４のいずれか１つに記載の方法を実行するために準備された機器。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、医療用画像の自動領域分割及び該自動領域分割用モデルの導出に関する。
【背景技術】
【０００２】
医療用の画像化は、診断及び患者の治療が改善される鍵として一般に認識されている。医療用の画像化はＸ線、コンピュータ断層撮影法（ＣＴ）、超音波及び磁気共鳴影像法（ＭＲＩ）等の画像診断法によって、最近数年間にわたって爆発的に発展してきた。従来、医療用画像は視覚的に検査され、重要な解剖学的構造又は一部の病巣（腫瘍）は、訓練を受けた放射線科医によって手作業でスケッチされている。しかしながら、このような処理は単調で時間がかかり、対象物の境界を正確にスケッチするために十分に注意して行われる必要がある。放射線科医の作業の一部又は全てを置き換えるためには、コンピュータ支援診断において、医療用画像データの領域分割が不可欠である。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】カワタ  ワイ  他（Kawata Y et al）、「曲率に基づく分析を用いた肺結節のコンピュータ支援鑑別診断（Computer aided differential diagnosis of pulmonary nodules using curvature based analysis）」、画像解析と画像処理１９９９会報（Image Analysis and Processing，1999. Proceedings）、アメリカ合衆国、電気電子技術者協会  コンピュータ社会（ＩＥＥＥ Computer Society）、４７０−４７５頁、１９９９年９月２７日
【発明の概要】
【発明が解決しようとする課題】
【０００４】
多くの場合、画像背景から対象物を分離することはとても難しいので、医療用画像の領域分割は厄介な仕事である。これは、処理における雑音が全ての医療用画像について固有である画像取得処理、及び、対象物自体の濃淡値マッピングの特性による。全ての画像取得装置の解像度には限度がある。それゆえ、医療用画像における各ボクセルの値は、幾つかの隣接する領域にわたって平均した値を示す。この現象は部分体積効果と呼ばれる。また、明暗差が小さい、大きさが小さい又は対象物が複雑な構造領域に位置する等の特徴を対象物が有している場合、自動領域分割することがより難しくなる。例えば、病巣（血管近傍の小結節、胸膜近傍の小結節又は結腸ポリープ等）の強度は、周囲組織（血管又は胸膜壁等）にとても類似している。この場合、従来の強度又はモデルに基づく方法では、対象物を適切に領域分割することができない。
【課題を解決するための手段】
【０００５】
本発明に係る３次元医療用画像を領域分割する方法は、医療用画像の特性の空間的な変動を示す医療用画像データを受け取るステップと、前記医療用画像データから形状データを導出するステップと、複数のモードを生成するために、前記医療用画像データ及び形状データを用いて前記医療用画像データをクラスタリングするステップと、前記複数のモードから導出され、前記医療用画像に関する領域分割データを出力するステップとを備えることを特徴とする。
【０００６】
本発明に係る３次元医療用画像を領域分割する方法は、医療用画像の特性の空間的な変動を示す医療用画像データを受け取るステップと、複数のモードを生成するために、前記医療用画像データをクラスタリングするステップと、前記複数のモードのうち、隣接するモードを結合するステップと、前記結合されたモードから導出され、前記医療用画像に関する領域分割データを出力するステップとを備えることを特徴とする。
【０００７】
本発明に係る３次元医療用画像を領域分割する方法は、医療用画像の複数のボクセルの強度を示す医療用画像データを受け取るステップと、該医療用画像データから前記複数のボクセル夫々に対する形状データを導出するステップと、前記複数のボクセルに対応し、夫々が対応するボクセルについて３次元座標、強度データ及び形状データを備える複数のベクトルを構成するステップと、複数のモードを生成するために前記複数のベクトルをクラスタリングするステップと、前記複数のモードのうち、隣接するモードを結合するステップと、前記結合した複数のモードから導出され、前記医療用画像に関する領域分割データを出力するステップとを備えることを特徴とする。
【０００８】
一の実施の形態では、領域分割方法は、自動的に医療の病巣を領域分割するために空間係数、強度係数及び体積形状係数をクラスタリングするステップを備える。提案するアルゴリズムは以下のステップ、（１）画像における各ボクセルに対する体積形状係数（ＳＩ）を計算するステップと、（２）５次元特徴ベクトルのセットを形成するためにＳＩ特徴を強度範囲及び空間位置（ｘ，ｙ，ｚ）と組み合わせるステップと、（３）５次元特徴ベクトルのセットをクラスタに分類するステップと、（４）隣接する領域又はモードを結合するために、クラスタ化アルゴリズムが提供する強度モードマップにおいて、空間的な特徴だけでなく形状的な特徴も考慮する改良期待値最大化アルゴリズム（ＥＭ）を用いるステップとを有する。結合した空間−強度−形状の特徴は、病巣又は腫瘍等の重要な解剖学的構造を領域分割するための優れた情報を提供する。
【０００９】
空間的な範囲（強度）情報だけでなく体積形状係数の特徴を考慮する５次元特徴ベクトルを生成する。即ち、領域分割された対象物における形状の特徴を計算する代わりに、領域分割の処理に形状の特徴を用いる。
【００１０】
隣接する領域又はモードを結合するために、クラスタ化アルゴリズムが提供する強度モードマップに改良期待値最大化アルゴリズムを適用する。事前確率の計算にフィルター処理された体積形状の特徴を提供するために形状係数モードマップを用いる。空間的な事前確率にこれらの形状特徴が考慮されるので、事前確率は空間的な情報だけでなく形状情報によっても決まる。結合した空間−強度−形状の特徴は重要な解剖学的構造又は病巣を領域分割するため優れた情報を提供する。重要な解剖学的構造又は病巣は、特に重要な領域が複雑な組織の領域内にあり、異なる形状の画像の強度が互いに非常に類似した（血管近傍の小結節、胸膜近傍の小結節又は結腸ポリープ等の）解剖学的構造又は病巣である。このような解剖学的構造又は病巣では、従来の強度又はモデルに基づく方法は対象物を適切に領域分割できない。
【００１１】
コンピュータで前述の方法を実施することが好ましい。本発明は、発明に係る方法を実行するための機器だけではなく、本発明に係る方法を実行するために用意されるコンピュータプログラム又はコンピュータプログラム製品にまで及ぶ。
【００１２】
本発明を、肺の小結節、肝臓病巣、マンモグラフィの腫瘤及び脳病巣等の様々な種類の病巣に適用することができる。更に、本発明の適用はＣＴスキャンに限定されず、ＭＲＩ画像、ＰＥＴ画像又はＸ線画像等の他の画像診断にも容易に適応することができる。
【図面の簡単な説明】
【００１３】
本発明の実施の形態を、下記に特定される図面を参照しながら説明する。
【図１】医療用画像装置と医療用画像装置からの画像データを処理するリモートコンピュータとを示す概要図である。
【図２】リモートコンピュータのより詳しい詳細図である。
【図３】本発明の実施の形態における方法の手順を示すフローチャートである。
【図４ａ】小結節及び血管を記録している２次元スライスの医療用画像である。
【図４ｂ】２次元スライスの医療用画像の形状係数マップである。
【図５ａ】試料の医療用画像である。
【図５ｂ】試料の医療用画像の強度モードマップである。
【図５ｃ】強度モードマップにおける小結線領域付近の強度モード値の図表である。
【図５ｄ】試料の医療用画像の形状係数モードマップである。
【図５ｅ】形状係数モードマップにおける小結線領域付近の形状係数モード値の図表である。
【発明を実施するための形態】
【００１４】
医療用の画像化
本発明はディジタル医療用画像に適用することができる。ディジタル医療用画像の一例はＣＴスキャン画像である。ＣＴスキャン画像は、ディジタル画像であり、人間又は動物の患者の部位をＣＴスキャンすることによって取得した一又は一連のＣＴ画像スライスから成る。各スライスは、スキャンされた部位におけるＸ線吸収の２次元ディジタルグレースケール画像である。スライスの特性は、用いられるＣＴスキャナによって決まる。例えば、高解像度のマルチスライス型ＣＴスキャナは、ｘ方向及びｙ方向（即ち、スライスの面）に０．５から１．０ｍｍ／画素の解像度で画像を生成することができる。各画素は３２ビットのグレースケールの解像度を有することができる。各画素の強度値をハウンズフィールド単位（ＨＵ）で表現することができる。連続するスライスをｚ方向（即ち、走査分離軸）に沿った一定の距離、例えば０．５から２．５ｍｍの距離で分離することができる。従って、スキャン画像は、スキャンされた部位及びスライス数によって決まる全体サイズを有する３次元（３Ｄ）グレースケール画像であってもよい。各画素は３次元空間でのボクセルであってもよい。また、スキャン画像は、単一のスライスであってもよいので、単一の２次元（２Ｄ）グレースケール画像であってもよい。
【００１５】
電子ビームコンピュータ断層撮像方法（ＥＢＣＴ）、マルチスライス法、スパイラルスキャン、又は、Ｘ線吸収を示す２Ｄ画像若しくは３Ｄ画像を出力として生成する技術等のＣＴスキャン技術によって、ＣＴスキャンを入手してもよい。
【００１６】
本発明はＣＴスキャン画像に限定されず、本発明をＭＲＩ画像、超音波画像又はＸ線画像等の他のディジタル医療用画像にも適用することができる。ディジタル化する前に、従来のＸ線画像をＸ線フィルムに現像してもよい。
【００１７】
コンピュータシステム
図１に示すように、スキャナ１０２からスキャン画像を受け取ってスキャン画像を構築するコンピュータ１０４はスキャン画像を生成することができる。スキャン画像は電子ファイル又は一連のファイルとして保存され、該電子ファイル又は一連のファイルは、固定式又は着脱可能式のディスク等の記憶媒体１０６に記憶される。スキャン画像はスキャン画像に関するメタデータを含んでもよい。コンピュータ１０４はスキャン画像を解析することができ、又は、例えば後述する、スキャン画像を処理するソフトウェアを実行するもう１つのコンピュータ１０８にスキャン画像を転送することができる。ソフトウェアは、着脱可能式ディスク若しくは固体メモリ等のキャリアに記憶されてもよく、又はローカルエリアネットワーク（ＬＡＮ）、広域エリアネットワーク（ＷＡＮ）若しくはインターネット等のネットワークを介してダウンロードされてもよい。
【００１８】
ここで述べたコンピュータは、図２に示すコンピュータシステム２００を構成することができる。本発明の実施の形態を、コンピュータシステム２００が実行するためのプログラム可能なコードとして実施することができる。本発明の様々な実施の形態を、本例のコンピュータシステム２００に関連して述べる。この説明を読めば、他のコンピュータシステム及び／又はコンピュータアーキテクチャを用いて本発明を実施する方法が当業者に明らかになる。
【００１９】
コンピュータシステム２００は、プロセッサ２０４等の一又は複数のプロセッサを有する。プロセッサ２０４は、特殊目的のディジタルシグナルプロセッサ又は汎用目的のディジタルシグナルプロセッサを含むが、これに限定されず、任意の種類のプロセッサであってもよい。プロセッサ２０４は通信基盤２０６（例えば、バス又はネットワーク）に接続される。様々なソフトウェアの実施を本例のコンピュータシステムに関連して述べる。この説明を読めば、他のコンピュータシステム及び／又はコンピュータアーキテクチャを用いて本発明を実施する方法が当業者に明らかになる。
【００２０】
コンピュータシステム２００は、主記憶装置２０８を更に有する。主記憶装置２０８はランダムアクセスメモリ（ＲＡＭ）であることが好ましい。コンピュータシステム２００は、更に補助記憶装置２１０を有してもよい。補助記憶装置２１０は、例えばハードディスクドライブ２１２、及び／又は、着脱可能な記憶ドライブ２１４を有することができる。着脱可能な記憶ドライブ２１４は、フレキシブルディスクドライブ、磁気テープドライブ又は光ディスクドライブ等に代表される。着脱可能な記憶ドライブ２１４は、公知の方法で着脱可能な記憶部２１８に対して読み書きする。着脱可能な記憶部２１８は、着脱可能な記憶ドライブ２１４によって読み書きされるフレキシブルディスク、磁気テープ又は光ディスク等に代表される。前述の内容からわかるように、着脱可能な記憶部２１８は、コンピュータが使用可能な記憶媒体を有する。コンピュータが使用可能な記憶媒体は、自身にソフトウェア及び／又はデータを記憶している。
【００２１】
別の実施では、補助記憶装置２１０はコンピュータシステム２００にコンピュータプログラム又は他の指示をロードさせる他の同様の手段を有してもよい。このような手段は、例えば、着脱可能な記憶部２２２及びインターフェース２２０を有する。このような手段の例として、（テレビゲーム装置で既に見られるような）プログラムカートリッジ及びカートリッジインターフェースと、（ＥＰＲＯＭ、ＰＲＯＭ又はフラッシュメモリ等の）着脱可能なメモリチップ及び関連するソケットと、他の着脱可能な記憶部２２２及び着脱可能な記憶部２２２からコンピュータシステム２００にソフトウェア及びデータを転送させるインターフェース２２０とを有してもよい。また、コンピュータシステム２００のプロセッサ２０４を用いて、着脱可能な記憶部２２２からプログラムの実行及び／又はデータのアクセスを行うことができる。
【００２２】
コンピュータシステム２００は、更に通信インターフェース２２４を有することができる。通信インターフェース２２４はソフトウェア及びデータをコンピュータシステム２００と外部装置との間で転送する。通信インターフェース２２４は、例として、モデム、（イーサネットカード（登録商標）等の）ネットワークカード、通信ポート、パーソナルコンピュータメモリカード国際協会（ＰＭＣＩＡ）スロット及びＰＭＣＩＡカード等を有してもよい。ソフトウェア及びデータは、通信インターフェース２２４を介して、信号２２８の形式で転送される。信号２２８は電子信号、電磁信号、光信号、又は、通信インターフェース２２４が受け取ることができる他の信号であってもよい。これらの信号２２８は、通信経路２２６を介して通信インターフェース２２４に与えられる。通信経路２２６は信号２２８を伝搬し、ワイヤ若しくはケーブル、光ファイバ、電話線、無線リンク、携帯電話リンク、高周波リンク、又は、任意の他の適切な通信チャンネルを用いて通信経路２２６を実装することができる。例えば、チャンネルを組み合わせて通信経路２２６を実装することができる。
【００２３】
このような適用では、「コンピュータプログラム媒体」及び「コンピュータで使用可能な媒体」の用語は、着脱可能な記憶ドライブ２１４、ハードディスクドライブ２１２に設置されたハードディスク及び信号２２８等の媒体に言及するために一般的に使用される。これらのコンピュータプログラム製品はソフトウェアをコンピュータシステム２００に提供する手段である。
【００２４】
（コンピュータ制御ロジックとも呼ばれる）コンピュータプログラムは、主記憶装置２０８及び／又は補助記憶装置２１０に記憶されている。また、コンピュータプログラムを、通信インターフェース２２４を介して受け取ることができる。このようなプログラムは、実行された場合、コンピュータシステム２００がここで論じる本発明を実施できるようにする。
【００２５】
従って、このようなコンピュータプログラムはコンピュータシステム２００の制御部を示す。ソフトウェアを用いて本発明を実施した場合、ソフトウェアは、幾つかの実施例を提供するために、コンピュータプログラム製品に記憶されて、着脱可能な記憶ドライブ２１４、ハードディスクドライブ２１２又は通信インターフェース２２４を用いてコンピュータシステム２００にロードされてもよい。
【００２６】
別の実施の形態では、本発明をハードウェア、ファームウェア、ソフトウェア、又は、これらの任意の組み合わせにおける制御ロジックとして実施することができる。
【００２７】
本発明の実施の形態に係る方法を説明する。特に定めのない限り、本発明の実施の形態に係る方法を、コンピュータシステム２００等のコンピュータシステムにおいて、ソフトウェアで実施すると理解される。「自動的」に実行されるステップへの言及は、好ましくは使用者の介在なしに、このようなソフトウェアによる実施を包含してもよい。
【００２８】
領域分割の要旨
図３は本発明の実施の形態に係る領域分割方法における手順を示している。
ステップＳ１にて、医療用画像において一組のボクセルの３次元位置を示すディジタルデータの形式で、各ボクセルの画像強度データと共に画像を入力する。
【００２９】
ステップＳ２にて、各ボクセルについて形状係数値を計算する。形状係数値は、各ボクセルで３次元形状を表す体積形状係数値であることが好ましく、各ボクセルでの強度の絶対値に無関係であってもよい。
【００３０】
ステップＳ３にて、各ボクセルの３次元座標を、各ボクセルの強度値及び形状係数値と組み合わせてボクセルに対応する一組の５次元ベクトルを生成する。
【００３１】
ステップＳ４にて、ベクトルをモードにクラスタリングする。これにより、モード内でのベクトルの類似性はモード間での類似性よりも高くなる。
【００３２】
ステップＳ５にて、例えば、空間的特徴及び形状特徴を考慮した改良期待値最大化法を用いて、モードを結合してセグメントにする。このステップでは、適切なレベルの領域分割を行うのに有益である。しかし、このステップは、ステップＳ４でのモード数が必要なセグメント数を満たすなら、状況に応じて選択的であってもよい。
【００３３】
ステップＳ６にて、導出した画像の領域分割結果を、例えば表示するものとして、又は、画像の更なる処理若しくは分析のための入力として出力する。
【００３４】
これらのステップ全てを、同じプロセッサ若しくはコンピュータシステムで、又は、同じソフトウェアアプリケーションによって実行する必要はない。例えば、ステップＳ１からＳ３の一又は複数のステップを、別々の前処理操作として実行し、その一方で、ステップＳ５からＳ６の一又は複数のステップを別々の後処理操作として実行することができる。
【００３５】
画像入力
ステップＳ１では、例えば、「医療用の画像化」及び／又は「コンピュータシステム」の節で前述したように、領域分割すべき画像又は必要な部分の画像を画像データとして入力する。各次元での概知のボクセル間隔から黙示的であるボクセルの座標を用いて、一組のボクセルに対応する強度値の３Ｄ行列として、画像データの形式を合わせる。また、対応する強度値と共に、セット中の各ボクセルの３次元座標を一組の４次元ベクトルとして画像データの形式を合わせてもよい。画像データを前処理して、例えば、雑音を削減したり、領域分割に不要な特徴を除いたりしてもよい。後述するように、ボクセルは、形状係数マップを計算させるために、十分に連続的又は隣接しているべきであるが、画像データは完全な一組の連続的なボクセルを示すことが本質ではない。
【００３６】
形状係数計算
ステップＳ２では、各ボクセルについて、ボクセル周辺の３Ｄの幾何学的特徴を示す形状係数値を計算する。形状係数値は、各ボクセルでの局所的な形状特徴を示す一方で、画像強度に無関係であることが好ましい。平面を除く全ての明確な形状が固有の形状係数に対応する。例えば、形状係数値１．００は球体の形状を示し、０．７５は円筒の形状を示す。体積形状係数は、同一表面を明示的に計算することなしに、各ボクセル周辺での同一表面の位相的な形状の特性を直接示す。
【００３７】
微分幾何学の第１及び第２基本形式に基づくガウス曲率及び平均曲率から形状係数値を計算することができる。各ボクセルｐで第１及び第２基本形式を計算する実際的なアプローチは、エイチ．ヨシダ（H. Yoshida）及びジェイ．ナッピ（J. Nappi）著、「結腸ポリープ検出のための３次元コンピュータ支援診断機構（Three-dimensional computer-aided diagnosis scheme for detection of colonic polyp）」、医療用画像化に関する電気電子技術者協会の報告書（IEEE Trans. Medical Imaging）、ｖｏｌ．２０、ｎｏ．１２、１２６１−１２７３頁、２００１年によって述べられるように、ｘ，ｙ，ｚの平滑化第１及び第２偏導関数を用いることである。
【００３８】
ある具体例では、各ボクセルでの体積形状係数値は以下のように定義される。
【００３９】
【数１】
2011526508.tif 000003
【００４０】
ここで、ｋ1 （ｐ）及びｋ2 （ｐ）は、ボクセルｐでの主曲率であり、以下のように定義される。
ｋ1 （ｐ）＝Ｈ（ｐ）＋√（Ｈ2 （ｐ）−Ｋ（ｐ）），
ｋ2 （ｐ）＝Ｈ（ｐ）−√（Ｈ2 （ｐ）−Ｋ（ｐ））      （２）
ここで、Ｋ（ｐ）及びＨ（ｐ）夫々はガウス曲率及び平均曲率である。
【００４１】
各ボクセルに関するガウス曲率及び平均曲率を、例えば、ジェイ．コーエンデリンク（J. Koenderink）及びエイ．ジェイ  バン  ドーン（A. J van Doorn）著、「表面形状と曲率スケール（Surface shape and curvature scales）」、画像及び映像の計算（Image and Vision Computing）、ｖｏｌ．１０、５５７−５６５頁、１９９２年に述べられているように、ヘッシアン行列から決定することができる。ヘッシアン行列は各ボクセル周辺での局所的な差分構造をとらえる。もう一つの代替手段では、ヘッシアン行列の固有値を直接用いて、形状係数値を計算する。
【００４２】
形状係数の特性を説明するために、図４ａ及び図４ｂ夫々は、２次元スライスの肺の原画像及び原画像に対応する形状係数マップを示している。以下の表１及び表２は、球状の小結節及び円筒状の血管夫々に対する形状係数値を示している。小結節に関する形状係数値の平均値は、血管に関する形状係数値の平均値よりも高い。
【００４３】
【表１】
2011526508.tif 000004
【００４４】
特徴ベクトルの形成
ステップＳ３にて、各ボクセルについて、形状係数値、画像強度値及びボクセルの３次元座標を組み合わせて、５次元空間Ｒ5 での５次元ベクトルを生成する。入力画像データと同様に、各ボクセルの３次元座標を、特徴ベクトルで明示的に示してもよく、又は３次元行列内のボクセルの位置から黙示的に示してもよい。例えば、各ボクセルにて、行列内でのボクセルの位置から導出される強度、形状係数及びボクセルの３次元座標によって５次元ベクトルを示すことができる。
【００４５】
クラスタリング
ステップＳ４にて、ベクトルをモードにクラスタリングする。クラスタリングは、特徴空間に基づくアルゴリズムを取り込み、データセットを幾つかのモード（別名クラスタ又はグループ）に分割する。これにより、モード内の類似性はモード間の類似性よりも大きくなる。特徴空間解析における画像領域分割の目的はこれらのモードの描写である。本発明の実施の形態に利用することができる別のノンパラメトリックなクラスタリング手法２つを以下の節、「平均シフトクラスタリング」及び「減法クラスタリング」で論じる。しかしながら、他のクラスタリング手法を用いてもよい。
【００４６】
平均シフトクラスタリング
平均シフト手法は特徴空間クラスタリングの非常に効果的な手法である。クラスタ数の予備的情報を必要とせず、クラスタの形状を制約しないのは、密度勾配のノンパラメトリック推定量である。平均シフト手法は、空間領域におけるボクセル位置等の空間情報と、強度レベル及び形状特徴等の特徴情報とを一緒に考慮する。平均シフト方法では、領域分割後に単に特徴を計算するというよりもむしろ、特徴を領域分割の処理に簡単に組み入れることができる。特徴誘導の領域分割を用いることによって、正確な対象物の境界スケッチを得ることができる。
【００４７】
ｄ次元空間Ｒd におけるｎ個のデータポイントｘi （ｉ＝１，．．．，ｎ）を考えた場合、カーネルＫ（ｘ）と窓半径ｈとを用いて得られる多変数カーネル密度推定量は以下のようになる。
【００４８】
【数２】
2011526508.tif 000005
【００４９】
放射対称性のカーネルにとって、以下の式を満たすカーネルｋ（ｘ）の外形を定義するのは十分である。
Ｋ（ｘ）＝ｃk,d ｋ（‖ｘ‖2 ）                          （４）
ここでｃk,d は、Ｋ（ｘ）の積分値が１になることを保証するように規格化された定数である。
【００５０】
密度関数のモードは、勾配関数のゼロ地点、即ち∇ｆ（ｘ）＝０で定義される。
密度推定量（３）式の勾配は以下の通りである。
【００５１】
【数３】
2011526508.tif 000006
【００５２】
ここでg(s)=-k’(s)である。第１項はカーネルＧ（ｘ）＝ｃg,d ｇ（‖ｘ‖2 ）を用いて計算したｘでの密度推定量に比例する。第２項は以下のようになる。
【００５３】
【数４】
2011526508.tif 000007
【００５４】
第２項は平均シフトである。平均シフトベクトルは常に密度の最大増加の方向を指す。
【００５５】
平均シフトアルゴリズムは以下の手順で構成される。
１）スタートポイントｘ0 を選択する。
２）平均シフトベクトルｍh （ｘt ）を計算する。
３）反復式ｘt+1 ＝ｘt ＋ｍh （ｘt ）に応じて、ｘt からｘt+1 に特徴ポイントを移動する。
４）‖ｘt −ｘt+1 ‖＜εであるなら、スタートポイントｘ0 からｘt までの経路上にある全てのポイントを、ｘt+1 を中心とするクラスタに分類する。
【００５６】
平均シフト手法は適合可能な勾配上昇手法の１つである。平均シフトベクトルの大きさは勾配そのものによって決まる。即ち、極大値に近い場合、ステップは小さく正確になる。結合した空間−強度の領域で作用することによって、平均シフトは、重要な強度変動を持つ可能性がある多くのボクセルをまとめて数個のモードにして、モードマップにする。
【００５７】
平均シフトクラスタリングにおける形状係数
５次元空間Ｒ5 におけるｎ個のデータポイントｘi （ｉ＝１，．．．，ｎ）（ここでｎは画像内のボクセル全ての数）を考えた場合、多変数カーネルは３つの放射対称性のカーネルの積で定義され、ユークリッド計量は以下に示す各カーネル関数について単一窓半径ｈを与える。
【００５８】
【数５】
2011526508.tif 000008
【００５９】
ここで、ｘs は空間位置、ｘr はグレーレベル範囲及びｘsiは形状係数特徴である。ｋ（ｘ）は全ての領域で用いられる共通の外形である。ｈs 、ｈr 及びｈsiは空間カーネル関数、範囲カーネル関数及び形状係数カーネル関数夫々に対するカーネル窓の大きさである。各カーネル窓の大きさを固定した定数に設定することができ、又は各カーネル窓の大きさを画像データに基づいて適合的に決めることができる。
【００６０】
平均シフト手法を用いることによって、形状係数の特徴を、対象物を領域分割するための強度特徴に組み入れることができる。この節で述べられている例は、平均シフトを用いた例に限定されない。即ち、減法クラスタリングを用いて、同様の結果を得ることができる。平均シフトクラスタリングは、特徴空間の基礎となる多変数配分のモードを推定する。クラスタ数は、空間で最も密度が濃い領域の中心（モード）を見つけることによって自動的に得られる。（６）式と同様に、３つのカーネル窓（空間カーネル窓、強度カーネル窓及び形状係数カーネル窓）を伴った平均シフトベクトルを以下のように計算することができる。
【００６１】
【数６】
2011526508.tif 000009
【００６２】
【数７】
2011526508.tif 000010
【００６３】
前述の手順は密度関数の勾配がゼロになるポイントｘit+1（モード）に収束するように保証されている。収束後、検出したモードはクラスタの中心である。
【００６４】
例として、図５ａは血管に付着した小結節の画像である。図５ｂ及び図５ｄ夫々は前述の結合した空間−形状係数−強度平均シフトの実施の形態を用いて生成された画像の強度モードマップ及び形状係数モードマップである。図５ｃ及び図５ｅは、図５ｂ及び図５ｄ夫々に対応する小結節領域周辺のモード値を示している。モード強度を、元のボクセルについて「フィルター処理した」値として確認することができ、モード強度値では異常値による品質劣化が少ない。
【００６５】
図５ｃに示すように、強度モードマップでモード値「−７２」によって示される小結節を背景（この場合、血管）から領域分割することができる。図５ｅでは、高い形状係数モード値（例えば８４）が小結節領域で得られる一方で、背景の値はかなり低い。
【００６６】
減法クラスタリング
減法クラスタリングは別のアプローチである。減法クラスタリングでは、データ点をクラスタの中心となる候補として考慮する。
【００６７】
ｄ次元空間Ｒd におけるｎ個のデータ点ｘi （ｉ＝１，．．．，ｎ）を考えた場合、各データ点がクラスタの中心となる候補であると仮定して、データ点ｘi での密度の大きさは以下のように定義される。
【００６８】
【数８】
2011526508.tif 000011
【００６９】
ここで、γa は正の定数であり、隣接領域を定義する。これにより、データ点が隣接するデータ点を多く有しているなら、データ点は高い密度値になる。
【００７０】
各データ点の密度の大きさを計算した後、最も高い密度の大きさを持つデータ点を最初のクラスタの中心として選択する。ｘc1が最初のクラスタの中心として選択されたデータ点である場合、Ｄc1はｘc1の密度の大きさである。次に、各データ点ｘi に関する密度の大きさは以下の式によって修正される。
【００７１】
【数９】
2011526508.tif 000012
【００７２】
ここで、γb は正の定数である。それゆえ、最初のクラスタの中心ｘc1近くのデータ点は、著しく減算された密度の大きさを有し、次のクラスタの中心として選択される可能性の低いデータ点となる。定数γb は、通常、クラスタの中心が密集することを防止するため、γa よりも大きい。一般的に、γb は１．５γa に等しい。
【００７３】
各データ点に関する密度の大きさを修正した後、次のクラスタｘc2を選択し、データ点に関する密度の大きさ全てを再び修正する。この処理を十分な数のクラスタの中心が生成されるまで繰り返す。
【００７４】
モード結合
クラスタリングアルゴリズムによって２つのモードマップ、即ち強度モードマップ（Ｍi ）と形状係数モードマップ（Ｍsi）とを生成する。各モードは特徴空間の任意の領域におけるデータの局所的構造を表す。モード数はカーネル窓の大きさ及びデータ構造によって決まる。この数は初期データから大きく圧縮されているけれども、この数は、まだ目的のクラス数よりも大きい。即ち、この場合には、過度の領域分割問題がある。過度の領域分割問題は、領域分割すべき目的の対象物を小さな領域に分割しすぎていることを意味する。
【００７５】
ステップＳ５で異なる／近いモードを結合するため、一の実施の形態では、グラフ理論手法が用いられている。グラフ理論手法は、ジェイ．シ（J. Shi）、ジェイ．マリク（J. Malik）著、「規格化された削減と画像領域分割（Normalized Cuts and Image Segmentation）」、パターン解析及び人工知能に関する電気電子技術者協会の報告書（IEEE Transactions on Pattern Analysis and Machine Intelligence）、ｖｏｌ．２２、ｎｏ．８、２０００年等に記載されている。このような手法を適合して、強度モードマップ及び形状係数モードマップから計算されたエッジの重みを含むように適合させることもできる。もう一つの実施の形態では、形状係数モードマップから空間情報及び形状情報を考慮し、強度モードマップＭi に用いることができる期待値最大化法（ＥＭ）を有する今までにない混合ガウスモデルが用いられている。ベイズ確率理論に基づいて、各モードについて、一つのクラスに属するモードの確率を以下のように定義する。
【００７６】
【数１０】
2011526508.tif 000013
【００７７】
ここで、ｐl （Ｍi ｜φl ）は、パラメータφl ＝（μl ，σl ）（平均値μl 及び標準偏差σl ）を用いたｌ番目のガウスモデルである。ｐ（φl(i)）は空間的な事前確率である。
【００７８】
各モードをガウス成分に割り当てる。ガウス成分について、各モードは以下の最も高い可能性を与える。
【００７９】
【数１１】
2011526508.tif 000014
【００８０】
空間的な事前確率ｐ（φl(i)）は（１１）式で重要な役割を果たす。以下に示すマルコフ確率場（ＭＲＦ）−ギブス確率場（ＭＲＦ−ＧＲＦ）によって空間的制限を課すことできる。
【００８１】
【数１２】
2011526508.tif 000015
【００８２】
ここでＮ（ｉ）はモードｉの隣接領域である。νc （ｌi ）は小集団ｃに関するポテンシャル関数であり、以下のように定義される。
【００８３】
【数１３】
2011526508.tif 000016
【００８４】
ここで、２次の隣接領域システム及び隣接領域間の対相互作用を考慮している。βは正の定数であり、クラスタリングの大きさを制御する。
【００８５】
強度モードマップと同様に、形状係数モードマップ（Ｍsi）は重要な形状情報を提供する。空間的な事前確率におけるこれらの形状特徴を考慮することによって、事前確率は空間情報だけではなく形状情報によっても決まる。空間形状に基づく事前確率を以下のように定義することができる。
【００８６】
【数１４】
2011526508.tif 000017
【００８７】
ここでＮc は全小集団の数であり、ｗsiは形状情報を考慮する重み因子である。
【００８８】
検討中のボクセルが隣接するボクセルと同様の形状である場合、隣接するボクセルと同じクラスに属する確率が高いボクセルに高い重みを割り当てると理解することができる。形状の特徴を（１５）式に組み入れることによって、事前確率は空間的な情報だけではなく、形状情報も考慮する。これは、強度モードマップ（Ｍi ）で空間的情報を考慮するだけの場合と比較してより良い領域分割を提供することができる。
【００８９】
領域分割の出力
ステップＳ６にて、ステップＳ４で識別したモード又はグループを出力して、ステップＳ１で入力された画像の領域分割を示すことができる。例えば、強度モードマップ及び／又は形状係数モードマップが画像内の対象物の領域分割を示すので、強度モードマップ及び／又は形状係数モードマップを出力してもよい。閾値をいずれかのモードマップに適用して、２値化領域分割マップを提供してもよいし、少なくとも示されたモードの数を減らしてもよい。別に又は加えて、ステップ６を実行する実施の形態において、ステップＳ５で得られる結合したモードを、ステップＳ６で出力してもよい。例えば、最も高い可能性をもたらすガウス成分を出力してもよい。
【００９０】
領域分割情報を表示用に出力して、例えば、画像で識別された様々な対象物を強調してもよいし、背景として識別された画像の部分を除去若しくは抑制してもよい。別に又は加えて、例えば、領域分割処理によって識別された対象物が病巣を示す可能性があるか否かを判定するために、領域分割結果を、画像を解析する処理に出力してもよい。
【００９１】
別の実施の形態
本発明を、肺小結節、肝臓病巣、マンモグラフィ腫瘤及び脳病巣等の様々な種類の病巣に適用することができる。所要の病巣のタイプに適した領域分割技術を用いてもよい。
【００９２】
本発明の適用はＣＴスキャンに限定されない。例えば、本発明の特徴をＭＲＩ画像、ＰＥＴ画像又はＸ線画像に適用してもよい。
【００９３】
本発明の別の実施の形態が、前述の説明を読むことで明らかになるかもしれない。それにもかかわらず、このような別の実施の形態は本発明の範囲に含まれる。
【図１】
2011526508.tif 000018
【図２】
2011526508.tif 000019
【図３】
2011526508.tif 000020
【図４ａ】
2011526508.tif 000021
【図４ｂ】
2011526508.tif 000022
【図５ａ】
2011526508.tif 000023
【図５ｂ】
2011526508.tif 000024
【図５ｃ】
2011526508.tif 000025
【図５ｄ】
2011526508.tif 000026
【図５ｅ】
2011526508.tif 000027
【国際調査報告】
2011526508.tif 000028
2011526508.tif 000029
2011526508.tif 000030