説明

医用画像処理装置および方法、並びにプログラム

【課題】3次元医用画像中の気管支のグラフ構造をより高い精度で構築する。
【解決手段】連続領域抽出部32が、気管支が表された3次元医用画像から、太い気管支内部の空気領域に対応する画素値を有する連続領域を抽出し、木構造構築部33が連続領域に対応する木構造を構築する。一方、線状構造抽出部34が、3次元医用画像中の各点の近傍の局所的な濃淡構造を解析することによって、細い気管支の断片を表す複数の線状構造を抽出する。木構造再構築部35は、太い気管支のグラフ構造を構成するノードと細い気管支の線状構造を表すノードを接続することによって、気管支全体を表すグラフ構造を再構築する。その際、太い気管支のノードと細い気管支の線状構造のノードとを結ぶ区間と、線状構造のノード同士を結ぶ区間とでは異なるコスト関数を用い、前者の区間におけるコスト関数は、画素値の変化が小さい区間ほど接続されやすくなるように定義した。

【発明の詳細な説明】
【技術分野】
【0001】
本発明は、3次元医用画像中の管状構造物を表すグラフ構造を構築する技術に関するものである。
【背景技術】
【0002】
肺癌は病気の進行により生存率が極端に低下していく病気の1つである。そのため、早期発見・早期治療が極めて重要である。例えば、胸部単純X線画像やCT画像などで腫瘍等と疑われる陰影が見られる場合、精密検査を行い、その陰影が本当に腫瘍なのか、本当に腫瘍である場合にはその良悪性について、判別を行う必要がある。この良悪性の判別には、気管支内視鏡を用いて腫瘍の一部を採取して病理検査を行う。この検査には、より迅速に、より正確に腫瘍の位置まで内視鏡を進めることが重要である。そのためには、検査前の撮影で得られたCT画像を用いて気管支の形状(分岐パターン等)と腫瘍までの経路とを把握することが有効である。
【0003】
ここで、CT等で得られた3次元医用画像から気管支等の線状構造を抽出する画像認識技術として、ヘッセ行列を用いた手法が提案されている。具体的には、まず、3次元医用画像に対して多重解像度変換を行った後、各解像度の画像に対してヘッセ行列の固有値解析を行い、線状構造要素を抽出する。この線状構造要素は、固有値解析によって得られる3つの固有値のうち1つだけが0に近いという特徴を有する。次に、各解像度の画像における解析結果を統合することによって、3次元医用画像中の様々なサイズの線状構造要素(血管)を抽出する。そして、最小全域木アルゴリズム等を用いて、抽出された各線状構造要素を連結することにより、3次元医用画像中の管状構造を表す木構造のデータが得られる。なお、最小全域木アルゴリズムによる線状構造要素の連結の際には、各線状構造要素の位置関係や、上記0に近い固有値に対応する固有ベクトルで表される、各線状構造要素の主軸方向に基づいたコスト関数が用いられる(特許文献1)。
【先行技術文献】
【特許文献】
【0004】
【特許文献1】特開2010−220742号公報
【発明の概要】
【発明が解決しようとする課題】
【0005】
一方、気管支は、気管から分岐を繰り返しながら、その直径が20mm以上のところから0.5mm以下のところまで徐々に細くなっていく木構造をしており、その太さによって異なる解剖学的・画像的特徴を呈する。具体的には、第1に分岐角度の違いが挙げられる。すなわち、気管支の太い箇所では分岐の角度が鈍角となりうるが、気管支の細い箇所では分岐の角度は鋭角となる。第2に、気管支の途切れが挙げられる。すなわち、気管支の太い箇所では疾患等による狭窄が生じていない限り、画像上で気管支が途切れてしまうことはないが、気管支の細い箇所では、パーシャルボリューム効果や心臓のモーションアーチファクトに起因する画像中のノイズにより、画像上では、気管支が途切れたように表されることがあり得る。第3に、気管支壁の厚さの違いが挙げられる。すなわち、気管支の太い箇所では、気管支壁が厚いため、気管支壁と気管支内部の空気領域とのCT値の差が顕著であるが、気管支の細い箇所では、気管支壁が薄いため、気管支壁と気管支内部の空気領域とのCT値の差が小さくなり、両者の境界が曖昧になる。
【0006】
上記特許文献1に記載の手法を用いて気管支を抽出する場合、各線状構造要素の主軸方向を加味したコスト関数を用いて最小全域木アルゴリズムを実行するため、気管支の細い箇所において、画像上の途切れがあったり、気管支壁と空気領域との画像上の差異が小さくなっていたりしても、各線状構造要素を正しく接続することが可能になる。逆に、気管支の太い箇所にある分岐の角度が鈍角の部分では、近接する線状構造要素間で主軸方向の差が大きくなるため、正しく接続されない可能性があり得る。
【0007】
本発明は上記事情に鑑みてなされたものであり、3次元医用画像中の気管支のグラフ構造をより高い精度で構築することができる医用画像処理装置および方法、並びにプログラムを提供することを目的とするものである。
【課題を解決するための手段】
【0008】
本発明の医用画像処理装置は、気管支が表された3次元医用画像を入力として、該気管支内部の空気領域に対応する画素値を有する連続領域を抽出する連続領域抽出手段と、該連続領域に対応するグラフ構造を構築するグラフ構造構築手段と、前記3次元医用画像を入力として、該3次元医用画像中の各点の近傍の局所的な濃淡構造を解析することによって、前記気管支の断片を表す複数の線状構造を抽出する線状構造抽出手段と、所定のコスト関数を用いて、前記グラフ構造を構成するノードと前記線状構造を表すノード、および、前記線状構造を表すノード同士を接続することによって、前記気管支を表すグラフ構造を再構築するグラフ構造再構築手段とを設け、前記グラフ構造を構成するノードと前記線状構造を表すノードとを結ぶ第1の区間と、前記線状構造を表すノード同士を結ぶ第2の区間とでは、異なるコスト関数を定義し、前記第1の区間におけるコスト関数を、該第1の区間のうち画素値の変化が小さい区間ほど接続されやすくなるように定義したことを特徴とする。
【0009】
本発明の医用画像処理方法は、気管支が表された3次元医用画像を入力として、該気管支内部の空気領域に対応する画素値を有する連続領域を抽出するステップと、該連続領域に対応するグラフ構造を構築するステップと、前記3次元医用画像を入力として、該3次元医用画像中の各点の近傍の局所的な濃淡構造を解析することによって、前記気管支の断片を表す複数の線状構造を抽出するステップと、所定のコスト関数を用いて、前記グラフ構造を構成するノードと前記線状構造を表すノード、および、前記線状構造を表すノード同士を接続することによって、前記気管支を表すグラフ構造を再構築するステップとを有し、前記グラフ構造を構成するノードと前記線状構造を表すノードとを結ぶ第1の区間と、前記線状構造を表すノード同士を結ぶ第2の区間とでは、異なるコスト関数が定義されたものであり、前記第1の区間におけるコスト関数は、該第1の区間のうち画素値の変化が小さい区間ほど接続されやすくなるように定義されたものであることを特徴とする。
【0010】
本発明の医用画像処理プログラムは、上記医用画像処理方法をコンピュータに実行させるためのものである。
【0011】
本発明において、線状構造を抽出する際に、線状構造の主軸方向を算出するようにし、第2の区間におけるコスト関数を、第2の区間の両端のノードにおける主軸方向の変化が小さい区間ほど接続されやすくなるように定義してもよい。
【0012】
また、3次元医用画像から異なる解像度の複数の画像を生成するようにし、生成された複数の画像のうちのより低解像度の画像を入力として連続領域を抽出するようにし、生成された複数の画像のうちのより高解像度の画像を入力として線状構造の抽出や主軸方向の算出を行うようにしてもよい。
【0013】
本発明において、気管支内部の空気領域に対応する画素値を有する連続領域を抽出する際には、所定の条件を満たす程度に画素値が近似する、近接した画素を連結することによって、連続領域を抽出するようにしてもよい。具体例としては、グラフカット法、領域拡張法、レベルセット法等を用いることができる。
【0014】
また、線状構造の抽出や主軸方向の算出には、ヘッセ行列の固有値解析を用いることができる。
【0015】
また、グラフ構造の構築には、最小全域木アルゴリズムを用いることができる。
【発明の効果】
【0016】
本発明によれば、気管支の太い箇所と細い箇所の各々の解剖学的・画像的特徴に適合した2つの抽出方法を組み合わせることによって、3次元画像中の気管支を表すグラフ構造を高い精度で生成することが可能になる。
【0017】
具体的には、第1に、3次元医用画像から気管支内部を表す連続領域が抽出される。前述のとおり、気管支の太い箇所は、途切れのない領域として画像中に表されるので、たとえ鈍角の分岐があっても、高い精度で気管支の連続領域が抽出される。したがって、この連続領域から生成されるグラフ構造も気管支を高い精度で表現したものとなる。
【0018】
第2に、3次元医用画像中の各点の近傍の局所的な濃淡構造を解析することによって、気管支の断片を表す複数の線状構造が抽出される。前述のとおり、気管支の細い箇所は画像中では途切れてしまうことがあり得るが、この第2の抽出手法では、気管支を断片的な線状構造の集合として取り扱うことができるので、高い精度での抽出が行われる。
【0019】
そして、気管支の連続領域を表すグラフ構造を構成するノードと気管支の断片を表す線状構造を表すノードとを結ぶ第1の区間と、気管支の断片を表す線状構造を表すノード同士を結ぶ第2の区間とでは、異なるコスト関数を用い、第1の区間におけるコスト関数を、第1の区間のうち画素値の変化が小さい区間ほど接続されやすくなるように定義し、このコスト関数を用いて気管支のグラフ構造が再構築される。したがって、上記第1の抽出方法で抽出された気管支の太い箇所と、上記第2の抽出方法で抽出された気管支の細い箇所とを適切に接続することが可能になる。
【0020】
以上により、3次元医用画像中の気管支のグラフ構造がより高い精度で作成される。
【図面の簡単な説明】
【0021】
【図1】本発明の実施の形態となる医用画像処理装置が導入された医用画像診断システムの概略構成図
【図2】本発明の実施形態における医用画像処理機能を実現する構成および処理の流れを模式的に示したブロック図
【図3】連続領域抽出部の詳細な構成を示したブロック図
【図4】木構造構築部の詳細な構成を示したブロック図
【図5】線状構造抽出部の詳細な構成を示したブロック図
【図6】肺野領域を模式的に表した図
【図7A】太い気管支を表す線状構造の集合を模式的に表した図
【図7B】肺野領域から太い気管支を表す線状構造の集合を除いた領域を模式的に表した図
【図7C】太い気管支を表す連続領域を模式的に表した図
【図8】太い気管支を表す連続領域を細線化した様子を模式的に表した図
【図9】細い気管支を表す線状構造の集合を模式的に表した図
【図10】太い気管支を表すノード点の集合と細い気管支を表すノード点の集合とを合成して模式的に表した図
【図11】ノード間でのコストの算出方法の第2の具体例を説明するための図
【図12】ノード間でのコストの算出方法の第3の具体例を説明するための図
【図13】太い気管支と細い気管支の境界付近を模式的に表した図
【発明を実施するための形態】
【0022】
以下、本発明の実施の形態となる医用画像処理装置が導入された医用画像診断システムについて説明する。
【0023】
図1は、この医用画像診断システムの概要を示すハードウェア構成図である。図に示すように、このシステムでは、モダリティ1と、画像保管サーバ2と、画像処理ワークステーション3とが、ネットワーク9を経由して通信可能な状態で接続されている。
【0024】
モダリティ1には、被検体の検査対象部位を撮影することにより、その部位を表す3次元医用画像の画像データを生成し、その画像データにDICOM(Digital Imaging and Communications in Medicine)規格で規定された付帯情報を付加して、画像情報として出力する装置が含まれる。本実施形態では、モダリティ1としてCTを用いて、被検体である人体の胸部を体軸方向にスキャンすることによって3次元画像データを生成する場合について説明する。
【0025】
画像保管サーバ2は、モダリティ1で取得された医用画像データや画像処理ワークステーション3での画像処理によって生成された医用画像の画像データを画像データベースに保存・管理するコンピュータであり、大容量外部記憶装置やデータベース管理用ソフトウェア(例えば、ORDB(Object Relational Database)管理ソフトウェア)を備えている。画像保管サーバ2は、画像処理ワークステーション3からの検索要求に応じて、画像データベースの検索を行い、抽出された画像データを画像処理ワークステーション3に送信する。
【0026】
画像処理ワークステーション3は、読影者からの要求に応じて、モダリティ1や画像保管サーバ2から取得した医用画像データに対して画像処理(画像解析を含む)を行い、生成された画像を表示するコンピュータであり、CPU,主記憶装置、補助記憶装置、入出力インターフェース、通信インターフェース、入力装置(マウス、キーボード等)、表示装置(ディスプレイモニタ)、データバス等の周知のハードウェア構成を備え、周知のオペレーティングシステムや、後述する医用画像処理プログラム等がインストールされたものである。この画像処理ワークステーション3で行われる医用画像処理は、CD−ROM等の記録媒体からインストールされた医用画像処理プログラムを実行することによって実現される。また、このプログラムは、インターネット等のネットワーク経由で接続されたサーバの記憶装置からダウンロードされた後にインストールされたものであってもよい。
【0027】
画像データの格納形式やネットワーク9経由での各装置間の通信は、DICOM等のプロトコルに基づいている。
【0028】
本実施形態では、画像処理ワークステーション3には、気管支や肺を対象とした手術や内視鏡検査のシミュレーションやナビゲーション機能が実装されている。この機能では、CTで得られた胸部の3次元医用画像から気管支を自動的に抽出した後、実際の内視鏡の挿入経路に沿って、気管支の仮想内視鏡画像を生成し、順次表示を行う。本発明の気管支を表すグラフ構造(以下では木構造を例とする)を構築する医用画像処理は、この気管支の自動抽出処理として実装されている。すなわち、画像処理ワークステーション3は、本発明の医用画像処理装置として機能するものである。
【0029】
図2は、画像処理ワークステーション3の機能のうち、本発明の実施形態となる気管支領域抽出処理に関連する部分の処理ブロックおよびデータや処理の流れを模式的に表した図である。本発明の医用画像処理もこの気管支領域抽出処理に含まれている。図に示すように、本実施形態における気管支領域抽出処理は、多重解像度画像生成部31、連続領域抽出部32、木構造構築部33、線状構造抽出部34、木構造再構築部35、気管支領域抽出部36によって実現される。また、3次元医用画像V、3次元医用画像Vの低解像度成分(以下、低解像度画像)VL、3次元医用画像Vの高解像度成分(以下、高解像度画像)VH、太い気管支の連続領域RL、太い気管支の木構造のノード情報NLn(nは各ノードを識別する添え字)、細い気管支の線状構造EHm(mは各線状構造を識別する添え字)、気管支全体の木構造T、気管支領域Rは、各々、上記各処理部によって、画像処理ワークステーション3の所定のメモリ領域に対して読み書きされるデータである。
【0030】
次に、本実施形態における気管支抽出処理の流れに沿って、上記各処理部の処理内容の概要について説明する。
【0031】
まず、多重解像度画像生成部31に対する入力データとなる3次元医用画像Vは、モダリティ1(CT)で撮影され、画像保管サーバ2に保管されたものであり、画像処理ワークステーション3からの検索要求に応じて、画像保管サーバ2によるデータベース検索によって抽出され、画像処理ワークステーション3に送信されたものである。多重解像度画像生成部31は、この3次元医用画像Vをメモリ領域から読み込んで多重解像度変換を行い、低解像度画像VL、高解像度画像VHを生成し、所定のメモリ領域に格納する。
【0032】
連続領域抽出部32は、低解像度画像VLをメモリ領域から読み込み、気管支内部の空気領域に対応する画素値を有する連続領域RLを抽出し、所定のメモリ領域に格納する。ここで、連続領域RLは、低解像度画像VLから抽出された気管支を表す領域であるから、後述の高解像度画像VLから抽出される気管支を表す領域よりも気管支の太い箇所を表すものとなる。
【0033】
木構造構築部33は、太い気管支を表す連続領域RLをメモリ領域から読み込み、太い気管支を表す木構造を構築し、そのノードの情報NLnを所定のメモリ領域に格納する。
【0034】
一方、線状構造抽出部34は、高解像度画像VHをメモリ領域から読み込み、高解像度画像VH中の各点の近傍の局所的な濃淡構造を解析することによって、細い気管支の断片を表す複数の線状構造EHmを抽出し、所定のメモリ領域に格納する。
【0035】
木構造再構築部35は、太い気管支を表す木構造の各ノードの情報NLnと、細い気管支を表す各線状構造EHmをメモリ領域から読み込み、後述の所定のコスト関数を用いて、太い気管支を表す木構造のノードNLnと細い気管支の線状構造を表すノードEHm、および、線状構造を表すノードEHm同士を接続することによって、気管支を表す木構造を再構築し、木構造データTを所定のメモリ領域に格納する。
【0036】
気管支領域抽出部36は、木構造データTをメモリ領域から読み込み、気管支領域全体Rを抽出し、所定のメモリ領域に格納する。
【0037】
次に、各処理部の詳細について説明する。
【0038】
多重解像度画像生成部31は、3次元医用画像Vに対して公知の多重解像度変換処理を行い、ガウシアンピラミッドを生成する。これによって生成される画像が低解像度画像VLと高解像度画像VHである。生成される画像は3以上であってもよく、その場合には、予め設定された解像度等の基準に基づいて、生成された各画像を低解像度画像群と高解像度画像群に分類すればよい。なお、多重解像度変換処理を行う前に、前処理として、人体の各成分に対応するCT値に基づく閾値処理やラベリング処理を行うことによって、被検体である人体の外側や体表面の軟部組織や骨等を表す連結領域を削除し、残された領域である肺野領域を抽出しておいてもよい。図6は、後続の処理対象となる肺野領域を模式的に表したものである。なお、肺野領域の抽出処理は、他の公知の手法を用いることもできる。
【0039】
連続領域抽出部32は、図3に示したように、ヘッセ行列固有値解析部32aとグラフカット処理部32bとからなる。
【0040】
まず、ヘッセ行列固有値解析部32aが、ヘッセ行列の固有値解析を行うことによって、気管支を表す線状構造を抽出する。図7Aは、抽出された気管支の線状構造の集合を模式的に表したものであり、図7Bは、図6の肺野領域から図7Aの気管支の線状構造の集合を除いた部分を模式的に表したものである。
【0041】
ここで、ヘッセ行列は、2階の偏微分係数を要素とする行列であり、3次元の画像に対しては、以下の3×3行列となる。
【数1】

【0042】
そして、以下のガウシアンカーネル(f)関数を用いた場合、ヘッセ行列を得るためのフィルタ係数は以下のように求められる。このとき、σは検出したいサイズの線状構造に対応させる。
【数2】

【0043】
このヘッセ行列の3つの固有値λ0、λ1、λ2は、画像中の線状構造を表す点では、以下に示すように、3つの固有値のうちの2つの固有値が大きく、残りの1つの固有値が0に近いという特徴を示す。
【数3】

【0044】
したがって、低解像度画像VL中の各点のうち、ヘッセ行列の3つの固有値が上記の特徴を有する点を線状構造として抽出することができる。
【0045】
なお、多重解像度画像生成部31が3以上の解像度の画像を生成する場合には、低解像度画像群に属する複数の低解像度画像の各々に対してヘッセ行列の固有値解析を行う。このとき、上記σを1種類のフィルタサイズに対応させた値としても、処理対象の画像の解像度が異なるため、複数のサイズの線状構造を抽出することが可能になる。そして、各解像度の画像における解析結果を統合することによって、最終的に、低解像度画像群中の様々なサイズの線構造を抽出することができる(Y Sato, et al.、「Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images.」、Medical Image Analysis、1998年6月、Vol.2、No.2、p.p.143-168等参照)。
【0046】
次に、グラフカット処理部32bは、公知のグラフカット法(例えば、特開2009-211138号公報等参照)を用いて、抽出された線状構造を表す点の集合(図7A参照)と、それ以外の点の集合(図7B参照)とを、前者の点が属する抽出対象領域と後者の点が属する背景領域とに分割することによって、前者の領域である太い気管支を表す連続領域RLを抽出する。図7Cは、抽出された太い気管支を表す連続領域RLを模式的に表したものである。
【0047】
なお、連続領域抽出部32は、気管支内部の空気領域に対応する画素値を有する連続した領域を抽出する手法(特に、所定の条件を満たす程度に画素値が近似する、近接した画素を連結する手法)であれば、領域拡張法やレベルセット法等の公知の他の手法を用いて、太い気管支を表す連続領域RLを抽出するようにしてもよい。
【0048】
木構造構築部33は、図4に示したように、細線化処理部33aと最小全域木処理部33bとからなる。
【0049】
まず、細線化処理部33aは、太い気管支を表す連続領域RLに対して、3次元の細線化処理を行う(詳細は、例えば、河田佳樹、外2名、「コーンビームCTの3次元血管像処理アルゴリズムについて」、電子情報通信学会論文誌、D-II Vol.J79-D-II No.6 pp.1134-1145 1996年6月等参照)。図8は、細線化された太い気管支を表す連続領域を模式的に表したものである。
【0050】
次に、最小全域木処理部33bは、連続領域を細線化することで得られた各画素をノード点として設定し、これらの点のうち、気管の位置に最も近い点を木構造の根(ルート)に決定する。この気管の位置に最も近い点は、公知の機械学習等によって得られた判別器を用いて決定するようにしてもよいし、細線化された連続領域のうち体軸方向において頭部に最も近い点に決定するようにしてもよいし、ユーザの手動操作によって指定するようにしてもよい。
【0051】
そして、最小全域木処理部33bは、複数のノード点の位置情報、すなわち各ノード点間の距離に基づいて各ノード点を接続する際のコストを算出し、最小全域木アルゴリズムを用いて、各ノード点を接続することによって、太い気管支を表す木構造データNLnを生成する。
【0052】
なお、木構造構築部33は、視線化された連続領域内を、決定された根から26近傍の探索を行うことによって木構造を構築するようにしてもよい。
【0053】
線状構造抽出部34は、図5に示したように、ヘッセ行列固有値解析部34aと主軸方向ベクトル特定部34bとからなる。
【0054】
まず、ヘッセ行列固有値解析部34aが、連続領域抽出部32のヘッセ行列固有値解析部32aと同様にして、ヘッセ行列の固有値解析を行うことによって、高解像度画像VH中の各点の近傍の局所的な濃淡構造を解析し、細い気管支を表す線状構造EHmを抽出する。図9は、抽出された線状構造EHmを模式的に表したものである。なお、ヘッセ行列固有値解析部34aで用いられるフィルタサイズを、ヘッセ行列固有値解析部32aで用いられたフィルタサイズと同じにしても、処理対象の画像の解像度が異なるため、高解像度画像VHからは、低解像度画像VLの場合よりも細い気管支を表す線状構造が抽出される。また、連続領域抽出部32と線状構造抽出部34の両方が、3次元医用画像Vにおける同じ領域から重複して気管支(連続領域、線状構造)を抽出してしまうことを回避するため、ヘッセ行列固有値解析部34aは、連続領域抽出部32で抽出された連続領域RL内からは、線状構造EHmを抽出しないようにしている。
【0055】
木構造再構築部35は、低解像度画像VLから構築された太い気管支を表す木構造データの各ノード点NLnと、高解像度画像VHから抽出された線状構造を表すノード点EHmとを用いて、最小全域木アルゴリズムを用いて、太い気管支を表すノード点NLnと細い気管支を表すノード点EHmとを接続することによって、気管支全体を表す木構造データTを再構築する。図10は、図8に示された太い気管支を表すノード点NLnと図9に示された細い気管支を表すノード点EHmと合成して模式的に表したものであり、図に示された各点を接続することによって、気管支全体を表す木構造データTが再構築される。
【0056】
まず、木構造再構築部35は、ノード点Piや、その前後のノード点Pi-1、Pi+1における複数の特徴のうちの1つ以上を用いて、ノード間のコストθi,i+1(i, i+1)を計算する。ここで、ノード点Pi-1、Pi、Pi+1は、太い気管支を表すノード点NLnと細い気管支を表すノード点EHmを区別せずに表したものである。以下、コストの算出方法の4つの具体例を挙げる。
【0057】
第1の具体例は、接続しようとする2つのノード点の位置関係を表すコストを算出する方法である。具体的には、以下の式によって、ノード点Piと点Pi+1の距離DPi,Pi+1からコストを計算する。この方法は、木構造構築部33の最小全域木処理部33bにおいて、太い気管支を表す木構造データNLnを構築する際にも用いられる。
【数4】

【0058】
第2の具体例は、各ノード点における気管支の進行方向の変化を表すコストを算出する方法である。図11に模式的に示したように、ノード点Piが有している気管支の進行方向ベクトルと直線Pii+1とのなす角度α、ノード点Pi+1が有している気管支の進行方向ベクトルと直線Pii+1とのなす角度βの和からコストを計算する。その計算は以下の式で表される。ここで、ノード点Pi、Pi+1が有している気管支の進行方向ベクトルとしては、ヘッセ行列固有値解析部32a、34aで算出された3つの固有値のうち0に最も近い固有値に対応する固有ベクトル、すなわち、ヘッセ行列固有値解析部32a、34aで抽出された線状構造の主軸方向ベクトルを用いることができる。この主軸方向ベクトルは、ヘッセ行列固有値解析部32a、34aにおいて、太い気管支を表す木構造のノード点NLn毎に、または、細い気管支の線状構造を表すノード点EHm毎に、その線状構造の点の位置情報等とともに関連づけて所定のメモリ領域に格納しておけばよい。
【数5】

【0059】
第3の具体例も、各ノード点における気管支の進行方向の変化を表すコストを算出する方法である。図12に模式的に示したように、ノード点Pi-1とノード点Piが既に繋がっている状況において、ノード点Pi-1からノード点Piに向かうベクトルと、ノード点Piからノード点Pi+1に向かうベクトルのなす角度γからコストを計算する。その計算は以下の式で表される。
【数6】

【0060】
第4の具体例では、ノード点Pとノード点Pi+1を結ぶ直線上における濃度値の変化を表すコストを算出する方法である。具体的には、以下の式によって、ノード点PiとPi+1との間の区間における濃度値の変化GPi,Pi+1からコストを計算する。
【数7】

【0061】
ここで、濃度値の変化GPi,Pi+1の具体例としては以下のものが挙げられる。
(1) ノード点PiとPi+1の間の区間からサンプリングされた複数の点における濃度値と、ノード点Piと点Pi+1の濃度を表す代表値(両点の濃度値の平均値等)との差の絶対値の最大値
(2) ノード点PiとPi+1の間の区間における濃度値の分散値
(3) ノード点PiとPi+1の間の区間のサンプリング点をソートし、ソートした点列内の濃度値の差分値
(4) ノード点PiとPi+1の間の区間のサンプリング点集合(x1, x2,・・・,xk,・・・,xK)中に外れ値が存在するか判別するために、次式のように、各サンプリング点における濃度値xkの、サンプリング点集合内の平均濃度値μからの偏差を、サンプリング点集合における不偏標準偏差σで割った検定統計量τkを用いたもの
【数8】

(5) 外れ値検定として用いられるスミルノフ・グラブス検定(正規分布を仮定する)やトンプソン検定を用いたもの
【0062】
次に、木構造再構築部35は、接続しようとする2つのノード点の位置関係(ノード点間の距離)、2つのノード点における気管支の進行方向ベクトルの変化、接続しようとする2つのノード点間の濃度値の変化の3つの要素に対して、接続しようとするノード点の種類に応じて異なる重みを付けて、2つのノード点を接続する際の接続コストを算出する。
【0063】
具体的には、低解像度画像VLから構築された太い気管支を表す木構造データの各ノード点NLn同士を接続しようとする際には、2つのノード点の位置関係の重みを大きくして接続コストを算出する。また、高解像度画像VHから構築された細い気管支を表すノード点EHm同士を接続しようとする際には、2つのノード点の進行方向ベクトルの変化の重みを大きくして接続コストを算出する(詳細は、上記特許文献1参照)。一方、低解像度画像VLから構築された太い気管支を表す木構造データの各ノード点NLnと、高解像度画像VHから構築された細い気管支を表すノード点EHmとを接続しようとする際には、2つのノード点間の画素値の変化の重みを大きくして接続コストを算出する。
【0064】
図13は、太い気管支を表す各ノード点NLnと、細い気管支を表すノード点EHmとの境界付近を模式的に表したものである。図の○印は前者のノード点、△印は後者のノード点を表しており、さらに実際の気管支壁が描かれている。この例の場合、ノード点NL100からNL104までの接続は、太い気管支を表すノード点同士の接続であるから、各ノード点の位置関係に大きな重みづけがなされた接続コストに基づいて、隣接するノード点が接続されている。そして、ノード点EH101からEH109までの接続は、細い気管支を表すノード点同士の接続であるから、ノード間での進行方向ベクトルの変化に大きな重みづけがなされた接続コストに基づいて、右下方向の進行方向を有するノード点EH101からEH104が接続される。ここで、ノード点EH101からEH104の接続順は、各ノードの位置関係のコストに基づいて、隣接するノード点が接続されている。また、同様にして、右方向に進行方向を有するノード点EH105からEH109も順に接続されている。これに対して、ノード点EH100は、他のノード点EH101, ・・・, EH109とは進行方向がまったく異なり、接続コストが所定の閾値よりも大きいため、これらの他のノード点EH101, ・・・, EH109とは接続されない。
【0065】
これらに対して、太い気管支を表すノード点NL104と細い気管支を表すノード点との接続では、接続しようとする2つのノード点間の濃度値の変化に大きな重みづけがなされた接続コストが計算される。したがって、太い気管支を表すノード点NL104と細い気管支を表すノード点EH101の間には気管支壁が存在しないので、各ノード点間の濃度値の変化が小さいため、両ノード点が接続される。同様に、太い気管支を表すノード点NL104と細い気管支を表すノード点EH105の間にも気管支壁が存在しないので、各ノード点間の濃度値の変化が小さいため、両ノード点が接続される。一方、太い気管支を表すノード点NL104と細い気管支を表すノード点EH100の間には気管支壁が存在するので、各ノード点間の濃度値の変化が大きくなるため、両ノード点は接続されない。なお、例えば、太い気管支を表すノード点NL104と細い気管支を表すノード点EH102が接続されないのは、ノード点NL104とノード点EH102との位置関係が、ノード点NL104とノード点EH101との位置関係よりも離れているため、前者の接続コストの法が後者の接続コストよりも大きくなってしまうからである。
【0066】
以上のようにして、太い気管支を表すノード点NL100からNL104、太い気管支を表すノード点NL104と細い気管支を表すノード点EH101、太い気管支を表すノード点NL104と細い気管支を表すノード点EH105、細い気管支を表すノード点EH101からEH104、細い気管支を表すノード点EH105からEH109が正しく接続される。一方、ノード点EH100は誤抽出点として扱われ、どのノード点とも接続されない。
【0067】
なお、木構造再構築部35は、上記のように、太い気管支を表すノード点NLnと細い気管支を表すノード点EHmのすべてについて、最小全域木アルゴリズムを用いて再接続する処理を行うようにしてもよいが、太い気管支を表すノード点NLnについては、木構造構築部33での接続結果をそのまま用いてもよい。これにより、結果的に重複する、太い気管支を表すノード点NLnの接続処理を省略することができ、処理効率が向上する。
【0068】
気管支領域抽出部36は、再構築された気管支全体を表すの木構造データTの各ノードにおいて、進行方向に直交する断面を求め、各断面において、前述のグラフカット法等の公知のセグメンテーション手法を用いて気管支の輪郭を認識し、その輪郭を表す情報を木構造データTの各ノードに関連づける。これにより、気管支領域Rが抽出される。
【0069】
以上のように、本発明の実施形態によれば、連続領域抽出部32、木構造構築部33が気管支の太い箇所の解剖学的・画像的特徴に適合した手法で太い気管支を抽出するとともに、線状構造抽出部34が気管支の細い箇所の解剖学的・画像的特徴に適合した手法で細い気管支を抽出する。そして、木構造再構築部35が、太い気管支を表す各ノード点同士を接続する際には、2つのノード点の位置関係の重みが大きくなるように定義されたコスト関数に基づいて接続コストを算出し、細い気管支を表すノード点同士を接続する際には、2つのノード点の進行方向ベクトルの変化の重みが大きくなるように定義されたコスト関数に基づいて接続コストを算出する。これにより、気管支の太い箇所と細い箇所の各々の解剖学的・画像的特徴に適合したコスト関数に基づいて算出された接続コストに基づいて、各ノード点が適切に接続される。さらに、太い気管支を表すノード点と、細い気管支を表すノード点とを接続する際には、2つのノード点間の濃度値の変化の重みを大きくして接続コストを算出する。これにより、これらのノード点間に気管支壁が存在するかどうかを考慮することができ、異なる手法で抽出された太い気管支と細い気管支とが適切に接続される。
【0070】
以上により、3次元医用画像中の気管支を表す木構造全体を高い精度で構築することが可能になる。
【0071】
上記実施形態はあくまでも例示であり、上記のすべての説明が本発明の技術的範囲を限定的に解釈するために利用されるべきものではない。また、上記の実施形態におけるシステム構成、ハードウェア構成、処理フロー、モジュール構成、ユーザインターフェースや具体的処理内容等に対して、本発明の趣旨から逸脱しない範囲で様々な改変を行ったものも、本発明の技術的範囲に含まれる。
【0072】
例えば、上記実施形態では、人体の気管支を抽出対象としているが、気管支と同様の解剖学的・画像的特徴を有する管状構造の抽出にも、本発明は適用可能である。
【0073】
また、上記実施形態では、本発明の気管支抽出手法が、気管支の仮想内視鏡画像を用いた、手術や検査のシミュレーションやナビゲーション機能の一部として実装されているが、慢性閉塞性肺疾患(COPD)等の肺や気管支の疾患を画像から診断するための気管支解析機能にも、本発明は実装可能である。この気管支解析機能では、CTで得られた胸部の3次元医用画像から気管支を自動的に抽出した後、ユーザによって指定された気管支の一部分について、CPR画像を生成したり、気管支の径や壁厚の計測を行ったりすることができる。
【0074】
さらに、上記実施形態では、多重解像度画像生成部31を設けて、低解像度画像から太い気管支を抽出することによって処理負荷の軽減、処理の高速化を図ったが、多重解像度画像生成部31を設けずに、もとの3次元医用画像から、太い気管支、細い気管支を抽出するようにしてもよい。この場合、例えば、連続領域抽出部32と線状構造抽出部34とで、ヘッセ行列の固有値解析の際のフィルタサイズを変えればよい。
【符号の説明】
【0075】
31 多重解像度画像生成部
32 連続領域抽出部
32a ヘッセ行列固有値解析部
32b グラフカット処理部
33 木構造構築部
33a 細線化処理部
33b 最小全域木処理部
34 線状構造抽出部
34a ヘッセ行列固有値解析部
34b 主軸方向ベクトル特定部
35 木構造再構築部
36 気管支領域抽出部

【特許請求の範囲】
【請求項1】
気管支が表された3次元医用画像を入力として、該気管支内部の空気領域に対応する画素値を有する連続領域を抽出する連続領域抽出手段と、
該連続領域に対応するグラフ構造を構築するグラフ構造構築手段と、
前記3次元医用画像を入力として、該3次元医用画像中の各点の近傍の局所的な濃淡構造を解析することによって、前記気管支の断片を表す複数の線状構造を抽出する線状構造抽出手段と、
所定のコスト関数を用いて、前記グラフ構造を構成するノードと前記線状構造を表すノード、および、前記線状構造を表すノード同士を接続することによって、前記気管支を表すグラフ構造を再構築するグラフ構造再構築手段とを備え、
前記グラフ構造を構成するノードと前記線状構造を表すノードとを結ぶ第1の区間と、前記線状構造を表すノード同士を結ぶ第2の区間とでは、異なるコスト関数が定義されたものであり、
前記第1の区間におけるコスト関数は、該第1の区間のうち画素値の変化が小さい区間ほど接続されやすくなるように定義されたものであることを特徴とする医用画像処理装置。
【請求項2】
前記線状構造抽出手段は、前記線状構造を抽出するとともに、該線状構造の主軸方向を算出するものであり、
前記第2の区間におけるコスト関数は、該第2の区間の両端のノードにおける前記主軸方向の変化が小さい区間ほど接続されやすくなるように定義されたものであることを特徴とする請求項1記載の医用画像処理装置。
【請求項3】
前記3次元医用画像から異なる解像度の複数の画像を生成する多重解像度画像生成手段をさらに備え、
前記連続領域抽出手段が、生成された複数の画像のうちのより低解像度の画像を入力として、前記連続領域を抽出するものであり、
前記線状構造抽出手段が、生成された複数の画像のうちのより高解像度の画像を入力として前記の処理を行うものであることを特徴とする請求項1または2記載の医用画像処理装置。
【請求項4】
前記連続領域抽出手段は、グラフカット法、領域拡張法、レベルセット法のいずれか1つを用いて前記連続領域を抽出するものであることを特徴とする請求項1から3のいずれか1項記載の医用画像処理装置。
【請求項5】
前記線状構造抽出手段は、前記3次元医用画像中の各点についてのヘッセ行列の固有値解析を行うものであることを特徴とする請求項1から4のいずれか1項記載の医用画像処理装置。
【請求項6】
前記グラフ構造構築手段および前記グラフ構造再構築手段は、全域木アルゴリズムを用いてグラフ構造を構築するものであることを特徴とする請求項1から5のいずれか1項記載の医用画像処理装置。
【請求項7】
気管支が表された3次元医用画像を入力として、該気管支内部の空気領域に対応する画素値を有する連続領域を抽出するステップと、
該連続領域に対応するグラフ構造を構築するステップと、
前記3次元医用画像を入力として、該3次元医用画像中の各点の近傍の局所的な濃淡構造を解析することによって、前記気管支の断片を表す複数の線状構造を抽出するステップと、
所定のコスト関数を用いて、前記グラフ構造を構成するノードと前記線状構造を表すノード、および、前記線状構造を表すノード同士を接続することによって、前記気管支を表すグラフ構造を再構築するステップとを有し、
前記グラフ構造を構成するノードと前記線状構造を表すノードとを結ぶ第1の区間と、前記線状構造を表すノード同士を結ぶ第2の区間とでは、異なるコスト関数が定義されたものであり、
前記第1の区間におけるコスト関数は、該第1の区間のうち画素値の変化が小さい区間ほど接続されやすくなるように定義されたものであることを特徴とする医用画像処理方法。
【請求項8】
コンピュータに、
気管支が表された3次元医用画像を入力として、該気管支内部の空気領域に対応する画素値を有する連続領域を抽出するステップと、
該連続領域に対応するグラフ構造を構築するステップと、
前記3次元医用画像を入力として、該3次元医用画像中の各点の近傍の局所的な濃淡構造を解析することによって、前記気管支の断片を表す複数の線状構造を抽出するステップと、
所定のコスト関数を用いて、前記グラフ構造を構成するノードと前記線状構造を表すノード、および、前記線状構造を表すノード同士を接続することによって、前記気管支を表すグラフ構造を再構築するステップとを実行させる医用画像処理プログラムであって、
前記グラフ構造を構成するノードと前記線状構造を表すノードとを結ぶ第1の区間と、前記線状構造を表すノード同士を結ぶ第2の区間とでは、異なるコスト関数が定義されたものであり、
前記第1の区間におけるコスト関数は、該第1の区間のうち画素値の変化が小さい区間ほど接続されやすくなるように定義されたものであることを特徴とする医用画像処理プログラム。

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図4】
image rotate

【図5】
image rotate

【図13】
image rotate

【図6】
image rotate

【図7A】
image rotate

【図7B】
image rotate

【図7C】
image rotate

【図8】
image rotate

【図9】
image rotate

【図10】
image rotate

【図11】
image rotate

【図12】
image rotate


【公開番号】特開2012−223315(P2012−223315A)
【公開日】平成24年11月15日(2012.11.15)
【国際特許分類】
【出願番号】特願2011−92723(P2011−92723)
【出願日】平成23年4月19日(2011.4.19)
【出願人】(306037311)富士フイルム株式会社 (25,513)
【Fターム(参考)】