(57)【特許請求の範囲】
【請求項１】
中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部と、
前記算出された評価行列に基づいて、前記画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出する評価部とを備えたことを特徴とする画像処理装置。
【請求項２】
中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部と、
前記算出された評価行列を固有値解析して固有値および固有ベクトルを算出し、算出された固有ベクトルの向きに基づいて線状構造の主軸方向または面状構造の法線方向を判別する方向判別部とを備えたことを特徴とする画像処理装置。
【請求項３】
前記中空球を表す関数は、下記式（３）により表されるものであることを特徴とする請求項１または２記載の画像処理装置。
【数１】
0005748693.tif 000013
ただし、ｘ、ｙ、ｚは実空間における３次元座標であり、Ｒは中空球の半径である。
【請求項４】
前記フィルタリング部は、複数サイズの前記中空球を表す関数について、該各中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列をそれぞれ算出するものであることを特徴とする請求項１から３のいずれか１項記載の画像処理装置。
【請求項５】
前記フィルタリング部は、前記中空球を表す関数の二階偏微分行列を用いたフィルタリングを、フーリエ空間において実施するものであることを特徴とする請求項１から４のいずれか１項記載の画像処理装置。
【請求項６】
前記画像が医用画像であり、前記中空の構造物が気管支であることを特徴とする請求項１から５のいずれか１項記載の画像処理装置。
【請求項７】
中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出し、
前記算出された評価行列に基づいて、前記画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出することを特徴とする画像処理方法。
【請求項８】
コンピュータを、
中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部と、
前記算出された評価行列に基づいて、前記画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出する評価部として機能させることを特徴とする画像処理プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像中の中空の構造物の構造を判別するための画像処理装置および方法並びに画像処理方法をコンピュータに実行させるためのプログラム、および、画像中に含まれる中空の構造物の主軸方向または法線方向を判別するための画像処理装置に関するものである。
【背景技術】
【０００２】
近年、医療機器（例えば多検出器型ＣＴ等）の進歩により質の高い３次元画像が画像診断に用いられる。３次元画像は多数の２次元画像から構成され情報量が多いため、医師が所望の観察部位を見つけ診断することに時間を要する場合がある。そこで、注目する臓器を抽出しＭＩＰ、ＶＲ、ＣＰＲ等の表示を行うことにより、臓器全体や病変の視認性を高め診断の効率化を図ることが行われている。
【０００３】
一方、医用画像中の血管および骨を抽出する手法として、ヘッセ行列を用いたヘッセ解析が提案されている（非特許文献１参照）。ヘッセ解析は、ガウシアンカーネル等の所定のフィルタについての二次微分カーネルを用いることにより算出した、二階の偏微分係数を要素とするヘッセ行列の固有値を解析することにより、画像中の局所構造が点、線および面のいずれであるかを判別するものである。ヘッセ解析を用いることにより血管は線状構造物として、骨は面状構造物として判別することが可能である。
【０００４】
しかしながら、非特許文献１の手法は、線状構造物の周辺に他の構造物（周辺構造物）が位置している場合に、この周辺構造物の一部を誤って線状構造物として判別する場合があった。非特許文献２の方法は、非特許文献１に提案されたフィルタを球形状の内側を１とし、外側を０とする、中実球を表す関数（中実球モデル関数）をたたみ込んだものに改良することにより、ヘッセ解析における二階の偏微分係数を計算する範囲を球表面に限定し、周辺構造物のフィルタリング結果に及ぼす影響を小さくすることができる。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】A. F. Frangi et. al., “Multiscale vessel enhancement filtering”, Proceedings of MICCAI, pp130-137, 1998.
【非特許文献２】M. Law et. al., “Three Dimensional Curvilinear Structure Detection Using Optimally Oriented Flux”, Proceedings of ECCV, pp368-382, 2008.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、非特許文献２の手法による中実球モデル関数が表す中実球モデルは、内側が何らかの充填物で満たされた構造的な特徴を有するものであるため、医用画像に含まれる気管支のような構造物の内側が空洞となる中空の構造物とは構造的な特徴が一致しない。このため、中空の構造物を含む画像に特許文献２に開示された中実球モデル関数をたたみ込んだフィルタを用いる手法を適用しても、モデル関数が表す形状と中空の構造物との構造的な特徴の違いに起因して、フィルタリング結果に誤差などが生じやすく、ヘッセ解析の結果の精度に限界がある。
【０００７】
このため、フィルタリングの精度を向上して、フィルタリングの結果として得られる評価行列に基づいて、より正確に中空の構造物の局所形状の判別（点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の判別）を行うことが求められている。また、同様に、フィルタリングの精度を向上して、フィルタリングの結果として得られる評価行列に基づいて、線状構造の軸方向または面状構造の法線方向をさらに高精度に判別することが求められている。
【０００８】
本発明は、上記事情に鑑みなされたものであり、画像に含まれる中空の構造物の点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の判別精度を向上することを第１の目的とする。また、画像に含まれる線状構造の軸方向または面状構造の法線方向を良好に判別することを第２の目的とする。
【課題を解決するための手段】
【０００９】
本発明による第１の態様の画像処理装置は、中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部と、算出された評価行列に基づいて、画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出する評価部とを備えたことを特徴とするものである。
【００１０】
本発明による第１の態様の画像処理方法は、中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出し、算出された評価行列に基づいて、画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出することを特徴とするものである。
【００１１】
本発明による第１の態様の画像処理プログラムは、コンピュータを、中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部と、算出された評価行列に基づいて、画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出する評価部として機能させることを特徴とするものである。
【００１２】
また、本発明の第２の態様の画像処理装置は、中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部と、算出された評価行列を固有値解析して固有値および固有ベクトルを算出し、算出された固有ベクトルの向きに基づいて線状構造の主軸方向または面上構造の法線方向を判別する方向判別部とを備えたものである。
【００１３】
上記「中空の構造物」とは、構造物の内部の少なくとも一部が空洞になっているものであれば何でもよい。例えば、中空の構造物は気管支などの管状（線状）構造物であってもよく、中空の粒状（点状）の構造物であってもよく、中空の面状の構造物であってもよい。また、中空の構造物は、２次元の構造物であっても３次元の構造物であってもよい。
【００１４】
上記本発明による第１または第２の態様の画像処理装置において、前記中空球を表す関数は、下記式（３）により表されるものであることが好ましい。
【数１】
0005748693.tif 000002
ただし、ｘ、ｙ、ｚは３次元空間の座標変数、ｒはその極座標表現における動径であり、Ｒは中空球の半径とする。
【００１５】
本発明による第１および第２の態様の画像処理装置において、前記フィルタリング部は、複数サイズの前記中空球を表す関数について、該各中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列をそれぞれ算出するものであることが好ましい。
【００１６】
また、本発明の第１または第２の態様の画像処理装置において、フィルタリング部は、中空球を表す関数の二階偏微分行列を用いたフィルタリングを、フーリエ空間において実施するものであることが好ましい。
【００１７】
また、本発明の第１または第２の態様の画像処理装置において、画像が医用画像であり、中空の構造物が気管支であることが好ましい。
【発明の効果】
【００１８】
本発明の第１の態様によれば、中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部と、算出された評価行列に基づいて、各画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出する評価部とを備えたため、判別対象物の構造的特徴を好適に表す中空球を表す関数により画像のフィルタリングを行うことができ、結果として評価行列から算出される評価値の精度を向上することができる。このため、評価値に基づいて画像に含まれる中空の構造物の判別精度を向上することができる。
【００１９】
本発明の第２の態様によれば、中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部と、評価行列を固有値解析して固有値および固有ベクトルを算出し、算出された固有ベクトルの向きに基づいて線状構造の主軸方向または面状構造の法線方向を判別する方向判別部とを備えたため、判別対象物の構造的特徴を好適に表す中空球を表す関数により画像のフィルタリングを行うことができ、結果として評価行列から算出される固有値および固有ベクトルの値の精度を向上することができる。このため、固有値および固有ベクトルに基づいて中空の構造物の各位置において、線状構造の軸方向または面状構造の法線方向を良好に判別することができる。
【図面の簡単な説明】
【００２０】
【図１】本発明の第１の実施形態による画像処理装置の構成を示す概略ブロック図
【図２】多重解像度変換を説明するための図
【図３Ａ】中空球モデル関数を説明するための図（その１）
【図３Ｂ】中空球モデル関数を説明するための図（その２）
【図４】線状構造の固有値を説明するための図
【図５】面状構造の固有値を説明するための図
【図６】本発明の第１の実施形態において行われる処理を示すフローチャート
【図７】本発明の第２の実施形態による画像処理装置の構成を示す概略ブロック図
【図８】本発明の第２の実施形態において行われる処理を示すフローチャート
【発明を実施するための形態】
【００２１】
以下、図面を参照して本発明の実施形態について説明する。図１は本発明の実施形態による画像処理装置の構成を示す概略ブロック図である。なお、図１のような画像処理装置１の構成は、補助記憶装置（不図示）に読み込まれたプログラムをコンピュータ（例えばパーソナルコンピュータ等）上で実行することにより実現される。また、このプログラムは、ＣＤ−ＲＯＭ等の情報記憶媒体に記憶され、もしくはインターネット等のネットワークを介して配布され、コンピュータにインストールされることになる。
【００２２】
画像処理装置１は、例えばＸ線ＣＴ装置２により撮像された複数の２次元画像を用いて、３次元画像Ｍ０を生成し、この３次元画像Ｍ０に含まれる線状構造物と面状構造物とを自動的に分割するものであって、画像取得部１０、検出領域設定部２０、判別部３０、分割部４０、表示部５０および入力部６０を備える。
【００２３】
画像取得部１０は、例えばＸ線ＣＴ装置２により撮像された複数のＣＴ画像（２次元画像）を取得し、複数の２次元画像から３次元画像Ｍ０を生成する。なお、画像取得部１０は、ＣＴ画像のみならず、いわゆるＭＲＩ画像、ＲＩ画像、ＰＥＴ画像、Ｘ線画像等の２次元画像を取得するものであってもよい。
【００２４】
検出領域設定部２０は、まず３次元画像Ｍ０のボクセルサイズを等方化する。例えば、３次元画像Ｍ０のボクセルサイズが、３次元画像Ｍ０におけるＸ，Ｙ，Ｚ方向のそれぞれにおいて、０．３ｍｍ、０．３ｍｍ、０．６ｍｍであった場合、（Ｘ，Ｙ，Ｚ）＝（０．５，０．５，０．５）（ｍｍ）に等方化する。このように、３次元画像Ｍ０のボクセルサイズを等方化することにより、後述する判別部３０において、Ｘ，Ｙ，Ｚ方向のそれぞれに同一サイズのカーネルを適用することができるため、演算を簡易に行うことが可能となる。
【００２５】
検出領域設定部２０は、等方化後に３次元画像Ｍ０を多重解像度変換して図２に示すように、解像度が異なる複数の３次元多重解像度画像Ｍｓｉ（ｉ＝０〜ｎ）（ガウシアンピラミッド）を生成する。なお、ｉ＝０は３次元画像Ｍ０と同一解像度、ｉ＝ｎは最低解像度を表す。３次元多重解像度画像Ｍｓｉのボクセルサイズは、解像度が高い順に、（Ｘ，Ｙ，Ｚ）＝（０．５，０．５，０．５）、（１．０，１．０，１．０）、（２．０，２．０，２．０）…となる。
【００２６】
判別部３０は、フィルタリング部３２および評価部３４を備える。フィルタリング部３２は、ヘッセ行列（評価行列）を用いたヘッセ解析を行うために、３次元多重解像度画像Ｍｓｉのそれぞれに、中空球モデル関数とガウシアンカーネルを用いたフィルタリングを行う。すなわち、解像度が異なる３次元多重解像度画像Ｍｓｉのそれぞれに対して、同一サイズのフィルタカーネルを畳み込む。また、中空球モデル関数の表す中空球の半径Ｒとガウスカーネルのサイズ（σ）は、予備的な解析等の知見に基づいて適切な値が設定される。なお、中空球の半径Ｒはガウスカーネルのサイズよりも少なくとも大きい値が設定される。
【００２７】
解像度が異なる３次元多重解像度画像Ｍｓｉのそれぞれに対して、同一サイズのフィルタカーネル（例えばR=2.0(voxel)、σ＝0.5(voxel)とする）を畳み込むことにより、３次元画像Ｍ０に対して、実質的に異なるサイズのフィルタカーネルを適用することとなるため、異なるサイズの点状構造物、線状構造物（例えば、気管支）および面状構造物（例えば皮質骨等の骨）を検出することができる。言い換えると、複数サイズの半径Ｒの中空球を表す関数について、該各中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列をそれぞれ算出している。
【００２８】
以下、本実施形態におけるヘッセ解析について説明する。ヘッセ解析に使用されるヘッセ行列は３次元画像に対しては下記の式（１）に示すように３×３の行列となる。
【数２】
0005748693.tif 000003
【００２９】
上記式（１）で示すヘッセ行列の各要素Ｉｘｘ、Ｉｘｙ、Ｉｙｙ、Ｉｙｚ、Ｉｚｚ、Ｉｚｘは、対象画像の画像データを中空球モデル関数ｆ（ｒ）およびガウシアンカーネル関数ｇ（ｒ）の２階偏微分を用いてフィルタリング（コンボリューション演算）して算出する。
【００３０】
本実施形態では、このフィルタリング処理をフーリエ空間で行う。まず対象画像の画像データをフーリエ変換する。次いで、このフーリエ変換された対象画像（FT(image)）に、それぞれフーリエ変換された中空球モデル関数、ガウシアンカーネル関数、および２階偏微分を積算する。そして、積算結果を逆フーリエ変換することによってフィルタリング結果を取得する。なお、この取得されたフィルタリング結果は、実空間でコンボリューション演算をして得られたフィルタリング結果と同じものとなる。中空球モデル関数のフーリエ変換をＦ（ν）、ガウシアンカーネル関数のフーリエ変換をＧ（ν）としたとき、以上の関係は式（２）のように表される。
【数３】
0005748693.tif 000004
【００３１】
ここで、本実施形態における中空球モデル関数ｆ（ｒ）は、式（３）に示すようにデルタ関数δ（ｒ−Ｒ）により定義される。なお、フーリエ空間においてフィルタリングを行うため式（３）をフーリエ変換した式（４）が式（２）において用いられる。
【数４】
0005748693.tif 000005
【数５】
0005748693.tif 000006
ただし、式（３）のｘ、ｙ、ｚは実空間における３次元座標であり、Ｒは中空球の半径である。また３次元フーリエ空間の変数（周波数）の極座標表現をν＝（νｘ２＋νｙ２＋νｚ２）１／２で表す。
【００３２】
また、式（２）におけるフーリエ変換されたガウシアンカーネル関数Ｇ(ν)は次式で表される。
【数６】
0005748693.tif 000007
なお、このガウシアンカーネル関数Ｇ（ν）は、非特許文献１および２と同様にフィルタリングの際に対象画像の画素値の微分範囲に幅を持たせるために用いる。
【００３３】
また、式（２）において、（２πνｘ）ｌ×（２πνｙ）ｍ×（２πνｚ）ｎの部分はフーリエ空間における微分処理に相当し、ｌ＋ｍ＋ｎ＝２（０＜ｌ、０＜ｍ、０＜ｎ）となるように、各微分方向に対応する係数ｌ、ｍ、ｎをそれぞれ式（２）に代入することにより、各要素Ｉｘｘ、Ｉｘｙ、Ｉｙｙ、Ｉｙｚ、Ｉｚｚ、Ｉｚｘに対応する値をそれぞれ算出できる。例えば、３次元多重解像度画像Ｍｓｉの処理の対象となる画素において、ヘッセ行列におけるＸ方向における二階偏微分値である要素Ｉｘｘは、式（２）にＸ方向にｌ＝２を代入し、Ｙ方向およびＺ方向にそれぞれｍ＝ｎ＝０を代入することにより算出できる。また、３次元多重解像度画像Ｍｓｉの処理の対象となる画素において、ヘッセ行列における要素Ｉｘｙは、式（２）にＸ方向およびＹ方向のそれぞれにｌ＝ｍ＝１を代入し、Ｚ方向にｎ＝０を代入することにより算出できる。
【００３４】
このようにして算出したヘッセ行列を固有値分解して固有値を得たとき、線状構造は図４に示すように、３つのうち２つの固有値の値が大きく、１つが０に近い特徴を持つことが知られている。例えば、式（１）の固有値は、線状構造からなる対象組織に対して、式（６）に示す関係を有する。なお、固有値は｜λ１｜＜｜λ２｜＜｜λ３｜とする。
【数７】
0005748693.tif 000008
【００３５】
また、面状構造は図５に示すように、３つのうち１つの固有値の値が大きく、２つが０に近い特徴を持つことが知られている。例えば、式（１）の固有値は、面状構造からなる対象組織に対して、式（７）のような関係を有する。
【数８】
0005748693.tif 000009
【００３６】
また、点状構造は、３つの固有値の全てが大きい特徴を持つことが知られている。例えば、式（１）の固有値は、点状構造からなる対象組織に対して、式（８）のような関係を有する。
【数９】
0005748693.tif 000010
【００３７】
したがって、固有値から線状構造らしさ、面状構造および点状構造らしさを判別することができ、判別結果を用いて、３次元画像Ｍ０において、線状構造物である気管支領域、および面状構造物である骨領域を分割することができる。
【００３８】
中空球モデル関数ｆ（ｒ）について、図３Ａ、３Ｂを用いて、以下説明する。図３Ａは、Ｘ方向における中空球モデル関数ｆ（ｒ）のレスポンスを実線で示し、中空球モデル関数ｆ（ｒ）のｘ方向の２階偏微分のレスポンスを破線で示している。図３Ｂの中央部分は、中空球モデル関数ｆ（ｒ）のレスポンスを、ｚ軸方向に等間隔に切断したｘｙ平面図により示している。また、図３Ｂの右部分は、中空球モデル関数ｆ（ｒ）のｘ方向の２階偏微分のレスポンスをｚ軸方向に等間隔に切断したｘｙ平面図により示している。なお、図３Ｂ中央部分、右部分とも各レスポンスを表す平面図をｚ座標値の大きい順に垂直方向上から並べて示している。図３Ｂのｘｙ平面図において、明度が大きい（白い）ほどフィルタ係数が正の大きい値となり、明度が小さい（黒い）ほどフィルタ係数が負の小さい値（フィルタ係数の絶対値の大きい値）となる。
【００３９】
中空の構造物に対して図３Ａに示すような中空球モデル関数を適用してフィルタリングを行った場合には、中空の構造物の輪郭部分が中空球モデルの球壁（半径ｒ＝−Ｒ，Ｒとなる部分）とそれぞれ一致した場合、式（２）により得られるレスポンスが好適に大きくなるため、非常に精度よく判別を行うことができる。このため、検出対象構造物を横断する線分、例えば気管支などの管状構造物の直径などに対して半径Ｒが一致するように、３次元多重解像度画像Ｍｓｉの解像度を設定することが好ましい。
【００４０】
評価部３４は、フィルタリング部３２が算出したヘッセ行列を固有値分解し、３つの固有値λ１，λ２，λ３（ただし、｜λ１｜≦｜λ２｜≦｜λ３｜とする）を算出する。そして、下記の式（１０）、（１１）に示すように、３次元多重解像度画像Ｍｓｉの各画素における線状構造らしさの評価値Ｌ０（Lineness）および面状構造らしさの評価値Ｐ０（Planeness）を算出する。なお、先述の通り、評価部３４は、ヘッセ行列の固有値λ１，λ２，λ３とそれぞれ対応する固有ベクトルｅ１，ｅ２，ｅ３に基づいて、点状構造、線状構造、面状構造を判別することが可能であるが、必ずしも点状構造、線状構造、面状構造の全てを評価する必要はなく、要求される仕様に応じて、点状構造、線状構造、面状構造の一部のみを判別してもよい。本実施形態では、医用画像から線状構造および面状構造を抽出することを目的としているため、線状構造および面状構造についてのみ評価値を算出している。
【数１０】
0005748693.tif 000011
なお、式（９）、（１０）におけるａ〜ｈは定数である。また、ＲA，ＲB，ＲCは下記の式（１１）〜（１４）により算出する。また、Ｓ2ndは二階偏微分値のパワーであり、下記の式（１４）により算出する。
【数１１】
0005748693.tif 000012
【００４１】
なお、本実施形態においては、解像度が異なる多重解像度３次元画像Ｍｓｉにおいて、それぞれ線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０が算出される。算出された評価値Ｌ０，Ｐ０は、元の３次元画像Ｍ０の対応する画素位置の評価値となるが、すべての多重解像度３次元画像Ｍｓｉの対応する画素位置において評価値が算出されることとなる。
【００４２】
分割部４０は、判別部３０が算出した線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０に基づいて、３次元画像Ｍ０の気管支領域および気管支以外の領域を領域分割する。具体的には、気管支領域を対象領域、気管支領域以外の領域を背景領域に設定し、３次元画像Ｍ０内の全画素位置において所定画素サイズの判別領域を設定し、Graph Cut領域分割方法を用いて、判別領域を対象領域と背景領域とに分割する。なお、Graph Cut領域分割方法は、「Yuri Y. Boykov, Marie-Pierre Jolly, “Interactive Graph Cuts for Optimal Boundary and Region Segmentation of Objects in N-D images”, Proceedings of “International Conference on Computer Vision”, Vancouver, Canada, July 2001 vol.I, p.105-112.」に記載されている。本実施形態では、分割部４０は、本発明者による過去出願である特開２０１１−２０６５３１号公報に記載されたものと同じ方法により、判別部３０が算出した線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０に基づいて、３次元画像Ｍ０の気管支領域および気管支以外の領域を領域分割する。
【００４３】
表示部５０は、２次元画像または３次元画像等を表示するモニタ、ＣＲＴ画面等である。本実施形態においては、表示部５０に、対象領域として分割された線状構造をボリュームレンダリング表示することにより線状構造や面状構造の全体を概観し、その連続性を可視化することができる。
【００４４】
入力部６０は、キーボード、マウス等である。
【００４５】
次いで、本実施形態において行われる処理について説明する。図６は本実施形態において行われる処理を示すフローチャートである。まず、画像取得部１０が、Ｘ線ＣＴ装置２により撮像された２次元画像から３次元画像Ｍ０を生成する（ステップＳＴ１）。次に、検出領域設定部２０が、３次元画像Ｍ０を等方化するとともに多重解像度変換して解像度が異なる複数の３次元多重解像度画像Ｍｓｉを生成する（ステップＳＴ２）。
【００４６】
次いで、判別部３０のフィルタリング部３２が、３次元多重解像度画像Ｍｓｉ毎に、中空球モデル関数およびガウシアンカーネル関数を用いたフィルタリングを行い、各画素位置におけるヘッセ行列を算出する（ステップＳＴ３）。次いで、評価部３４が、算出された各ヘッセ行列に基づいて線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０を算出する（ステップＳＴ４）。
【００４７】
そして、分割部４０が上述したGraph Cut領域分割方法を用いて、３次元画像Ｍ０を対象領域（気管支領域）と背景領域とに分割する（ステップＳＴ５）。そして、表示部５０が分割された対象領域および背景領域をボリュームレンダリング表示し（ステップＳＴ６）、処理を終了する。
【００４８】
このように、本実施形態によれば、中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部３２と、算出された評価行列に基づいて、各画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出する評価部３４とを備えたため、判別対象物の構造的特徴を好適に表す中空球を表す関数により画像のフィルタリングを行うことができる。このため、好適なフィルタリングの結果得られた評価行列から算出される評価値に基づいて、画像に含まれる中空の構造物の点状構造、線状構造および面状構造のうち評価されたものの評価精度を向上することができる。
【００４９】
また、上記第１の実施形態においては、複数サイズの半径Ｒの中空球を表す関数について、該各中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列をそれぞれ算出し、この複数の評価行列に基づいて、それぞれ各画素位置における線状構造らしさおよび面状構造らしさの評価値を算出している。このため、各カーネルサイズに対応するサイズの構造物を好適に評価することができる。また、この評価値に基づいて、上記の第１の実施形態のように対象領域を分割した場合には、好適に複数サイズの中空の構造物を対象領域として分割することができる。
【００５０】
なお、異なるサイズの構造物に対してヘッセ解析を行う手法として、本実施形態のように複数の解像度の画像に対して、一定のカーネルサイズσおよび一定の中空球モデル関数の半径Ｒを適用してヘッセ解析を行う手法を適用してもよく、一定の解像度の画像に対して、カーネルサイズσおよび中空球モデル関数の半径Ｒを異ならせてヘッセ解析を行う手法を適用しても同効果が得られる。例えば、中空の構造物の長軸や短軸の長さなど、検出対象構造物を横断する線分の長さを試験的に得られたデータなどから予め取得しておき、検出対象構造物を横断する線分の長さと半径Ｒの大きさが一致するように、３次元多重解像度画像Ｍｓｉの解像度、あるいは、カーネルサイズσおよび半径Ｒを設定することが考えられる。
【００５１】
なお、計算負荷の抑制と計算処理の高速化のためには本実施形態のようにフーリエ空間においてフィルタリングを行うことが好適であるが、実空間において対象画像に中空球モデル関数ｆ（ｒ）の２階偏微分およびガウシアンカーネルを畳み込み積分（コンボリューション）することによりフィルタリング行ってもよい。
【００５２】
また、上記第１の実施形態においては、線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０の双方を算出しているが、線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０のいずれか一方のみを算出するようにしてもよい。
【００５３】
また、上記第１の実施形態においては、ヘッセ行列を固有値分解することにより取得した固有値を用いて線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０を算出しているが、固有値を得ることなく、ヘッセ行列の要素をそのまま用いて線状構造らしさの評価値および面状構造らしさの評価値を算出するようにしてもよい。
【００５４】
本発明の第２の実施形態として、方向判別のための画像処理装置を以下に説明する。図７は方向判別のための画像処理装置の機能ブロック図を表す。第２の実施形態における画像処理装置は、第１の実施形態における分割部４０と評価部３４と、表示部６０を備えておらず、代わりに方向判別部３５を備えている。この方向判別部３５による方向判別処理以外の処理である、画像取得部１０による画像取得処理と、検出領域設定部２０による検出領域設定処理と、フィルタリング部３２によるフィルタリング処理については、第１の実施形態と同じであり、各機能ブロックの機能も共通している。以下、第１の実施形態と異なる点を中心に説明し、第１の実施形態と同じ点については説明を省略する。
【００５５】
方向判別部３５は、フィルタリング部３２が算出した評価行列を固有値解析して固有値および固有ベクトルｅ１，ｅ２，ｅ３を算出し、算出された固有ベクトルの向きに基づいて線状構造物の主軸方向または面状構造の法線方向を判別する。詳細には、フィルタリング部３２が得たヘッセ行列の固有値λ１，λ２，λ３および固有値ベクトルｅ１，ｅ２，ｅ３を取得し、式（６）に示す関係を満たしていれば、固有値λ１に対応する固有ベクトルｅ１の指す方向を線状構造の軸方向と判別し、式（７）に示す関係を満たしていれば、固有値λ３に対応する固有ベクトルｅ３の指す方向を面状構造の法線方向と判別する。
【００５６】
図８は本実施形態において行われる処理を示すフローチャートである。まず、画像取得部１０が、Ｘ線ＣＴ装置２により撮像された２次元画像から３次元画像Ｍ０を生成する（ステップＳＴ１１）。次に、検出領域設定部２０が、３次元画像Ｍ０を等方化するとともに多重解像度変換して解像度が異なる複数の３次元多重解像度画像Ｍｓｉを生成する（ステップＳＴ１２）。
【００５７】
次いで、判別部３０のフィルタリング部３２が、３次元多重解像度画像Ｍｓｉ毎に、中空球モデル関数およびガウシアンカーネル関数を用いたフィルタリングを行い、各画素位置におけるヘッセ行列を算出する（ステップＳＴ１３）。次いで、方向判別部３５が、ヘッセ解析の結果に基づいて、各画素の局所的な固有値および固有値ベクトルを取得し、先述のように式（４）および式（５）に基づいて線状構造物の軸方向または面状構造物の法線方向を算出し（ステップＳＴ１４）、処理を終了する。
【００５８】
以上のように、第２の実施形態においても、第１の実施形態と同様に判別対象物の構造的特徴を好適に表す中空球を表す関数により画像のフィルタリングを行うことができ、結果として評価行列から算出される固有値および固有ベクトルの値の精度を向上することができる。このため、固有値および固有ベクトルに基づいて中空の構造物の各位置において、線状構造の軸方向または面状構造の法線方向を良好に判別することができる。上記第２の実施形態における画像処理装置により得られた線状構造の軸方向や面状構造の法線方向は、例えば気管支や大腸など管状構造のＣＰＲ(Curved Planer Reformation / Reconstruction)法によって３次元医用画像から再構成されたＣＰＲ画像を作成する際の軸方向として用いることができ、線状構造の軸方向や面状構造の法線方向を必要とする周知の種々の処理に好適に用いることができる。なお、本発明により得られた線状構造の軸方向や面状構造の法線方向は、例えば、「Armin Kanitsar, Dominik Fleischmann, Rainer Wegenkittl, Petr Felkel, Eduard Groller: CPR - Curved Planar Reformation. IEEE Visualization 2002.」など周知の種々のＣＰＲ画像生成手法に用いることができる。
【００５９】
また、上記第２の実施形態においても、複数サイズの半径Ｒの中空球を表す関数について、該各中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列をそれぞれ算出し、この複数の評価行列に基づいて、それぞれ各画素位置における固有値および固有ベクトルを算出している。このため、各カーネルサイズに対応するサイズの構造物について線状構造の軸方向や面状構造の法線方向を好適に判別することができる。
【００６０】
なお、方向判別部３５は、ヘッセ行列の固有値と固有ベクトルに基づいて、式（６）に基づく線状構造の軸方向、または、式（７）に基づく面状構造の法線方向のいずれか一方のみを判別してもよい。
【００６１】
また、上記各実施形態においては、画像が医用画像であり、中空の構造物が気管支であるため、フィルタリングに用いる中空球を表す関数と画像に含まれる中空の構造物の構造的な特徴が好適に一致するため、第１の実施形態においては非常に精度よく点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出することができる。また、第２の実施形態においては非常に精度よく線状構造の軸方向または面状構造の法線方向を判別することができる。
【００６２】
なお、上記各実施形態では線状構造物として気管支の判別を例に挙げたが、中空構造であれば他の構造物の判別にも適用できる。また、中空の構造物と、血管などの中実の構造物の両方を含む画像や、中空の構造物と、骨、皮膚、葉間胸膜等の面状構造物の両方を含む画像にも利用することができる。
【００６３】
なお、上記各実施形態においては、３次元画像Ｍ０に含まれる線状構造物および面状構造物を対象に各評価処理または方向判別処理を行っているが、２次元画像に含まれる線状構造物および面状構造物を対象に各評価処理または方向判別処理を行ってもよい。
【００６４】
また、上記各実施形態においては、Graph Cut領域分割方法を用いて３次元画像Ｍ０に含まれる線状構造物と面状構造物とを分割しているが、Watershedアルゴリズム等の他の領域分割手法を用いてもよいことはもちろんである。ここで、Watershedアルゴリズムは、画像の画素値情報を高度とみなした地形において水を満たしていく際に、異なるくぼみに溜まる水の間で境界ができるように画像を分割する手法である。このため、３次元画像Ｍ０上の評価値Ｌ０，Ｐ０に対し、適当な平滑化を行った後でWatershedアルゴリズムを実行することによって、線状構造物および面状構造物を分割することが可能である。
【００６５】
上記の各実施形態はあくまでも例示であり、上記のすべての説明が本発明の技術的範囲を限定的に解釈するために利用されるべきものではない。
【００６６】
この他、上記の実施形態におけるシステム構成、ハードウェア構成、処理フロー、モジュール構成、ユーザインターフェースや具体的処理内容等に対して、本発明の趣旨から逸脱しない範囲で様々な改変を行ったものも、本発明の技術的範囲に含まれる。
【符号の説明】
【００６７】
１      画像処理装置
２    Ｘ線ＣＴ装置
１０    画像取得部
２０    検出領域設定部
３０    判別部
３２    フィルタリング部
３４    評価部
３５    方向判別部
４０    分割部
５０    表示部
６０    入力部
【図１】
0005748693.tif 000014
【図２】
0005748693.tif 000015
【図３Ａ】
0005748693.tif 000016
【図４】
0005748693.tif 000017
【図５】
0005748693.tif 000018
【図６】
0005748693.tif 000019
【図７】
0005748693.tif 000020
【図８】
0005748693.tif 000021
【図３Ｂ】
0005748693.tif 000022