JP5833958B2 - 画像処理装置および方法並びにプログラム - Google Patents

Description

本発明は、画像中の構造物の構造を判別するための画像処理装置および方法並びに画像処理方法をコンピュータに実行させるためのプログラムに関するものである。

近年、医療機器（例えば多検出器型ＣＴ等）の進歩により質の高い３次元画像が画像診断に用いられる。３次元画像は多数の２次元画像から構成され情報量が多いため、医師が所望の観察部位を見つけ診断することに時間を要する場合がある。そこで、注目する臓器を抽出しＭＩＰ、ＶＲ、ＣＰＲ等の表示を行うことにより、臓器全体や病変の視認性を高め診断の効率化を図ることが行われている。

一方、医用画像中の血管および骨を抽出する手法として、ヘッセ行列を用いたヘッセ解析が提案されている（非特許文献１参照）。ヘッセ解析は、ガウシアンカーネル等の所定のフィルタについての二次微分カーネルを用いることにより算出した、２階の偏微分係数を要素とするヘッセ行列の固有値を解析することにより、画像中の局所構造が点、線および面のいずれであるかを判別するものである。ヘッセ解析を用いることにより血管は線状構造物として、骨は面状構造物として判別することが可能である。

しかしながら、非特許文献１の手法は、線状構造物の周辺に他の構造物（周辺構造物）が位置している場合に、この周辺構造物の一部を誤って線状構造物として判別する場合があった。非特許文献２の方法は、非特許文献１に提案されたフィルタを、球（内側を１とし、外側を０とする中実球）を表す関数（中実球モデル関数）の２階偏微分をたたみ込んだものに改良することにより、各画素のフィルタリング範囲を画素を中心とする球表面に限定し、周辺構造物のフィルタリング結果に及ぼす影響を小さくすることができる。

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.

しかしながら、非特許文献２の手法を用いて、様々な太さの血管を含む医用画像から線状構造である血管を判別すると、血管の内側に実際の血管よりも細い血管が誤って認識されてしまう場合がある。言い換えると、非特許文献２の手法で線状構造物を判別する場合に、構造物の輪郭部分のうち、中実球モデル関数の中実球の曲率半径より大きい曲率半径を有する輪郭部分で、中実球モデル関数の表す中実球の直径と略同一の直径を有する線状の構造物が存在すると誤って判別されてしまうという問題がある。

上記問題は、中実球モデル関数の中実球よりも曲率半径が大きい構造物の輪郭部分であればどこでも起こりうるため、本発明は、この問題に鑑みて、ヘッセ解析において中実球モデル関数を用いたフィルタリングを行う画像処理方法において、中実球モデル関数の表す中実球よりも曲率半径が大きい構造物の輪郭部分で生じる構造物の誤判別を防止することを目的とする。

本第１の発明による画像処理装置は、画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、フィルタリング部が、画像中の各画素位置において、中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、取得した１階偏微分ベクトルを固有ベクトルの方向に射影した値を用いて、中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行う補正部を備えたものであることを特徴とするものである。

本第１の発明による画像処理方法は、画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング工程と、算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価工程とを備えた画像処理方法であって、フィルタリング工程が、画像中の各画素位置において、中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、取得した１階偏微分ベクトルを固有ベクトルの方向に射影した値を用いて、中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行うものである。

なお、本第１の発明による画像処理方法を、コンピュータに実行させるためのプログラムとして提供してもよい。

上記、「中実球の半径と同じ半径を有する中空球」とは、中実球と中空球の半径が厳密に一致する場合だけでなく、中空球を表す関数の１階偏微分ベクトルを用いて、中実球を表す関数の各方向の２階偏微分の中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す効果が得られる範囲であれば、中空球の半径が中実球の半径より大きい場合や小さい場合も含む。上記「一方の位置の応答波形を打ち消す効果」を良好に得るために、中空球の半径と中実球の半径ができるだけ等しいことが好ましく、例えば、中空球の半径と中実球の半径の差が２割以下であることが好ましく、中空球の半径と中実球の半径の差が１割以下であることが好ましい。

本第１の発明による画像処理装置において、ヘッセ行列の固有値をλ₁、λ_２、λ_３とし、補正後の固有値をλ₁ ^’、λ_２ ^’、λ_３ ^’ とし、固有ベクトルをｅ_１＝（ｘ_１,ｙ_１ ,ｚ_１），ｅ_２＝（ｘ_２,ｙ_２ ,ｚ_２），ｅ_３＝（ｘ_３,ｙ_３ ,ｚ_３）とし、第１の半径を有する中空球を表す関数の１階偏微分ベクトルを（ρ_１,ρ_２,ρ_３）とすると、補正部は、下記式（１２）で算出されるρ_１ ^’,ρ_２ ^’,ρ_３ ^’および所定の係数αを用いて、下記式（１３）に示すように、固有値を補正することにより、一方の応答波形を打ち消す補正を行うものであることが好ましい。

本第１の発明による画像処理装置において、中空球を表す関数は、下記式（１０）により表されるものであることが好ましい。
ここで、ｘ、ｙ、ｚは３次元空間の座標、ｒはその極座標表現であり、Ｒ_４は中空球の半径とする。

本第２の発明による画像処理装置は、画像中の各画素位置において、第１の中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、フィルタリング部が、画像中の各画素位置において、第１の中実球の半径である第１の半径より大きい第２の半径を有する第２の中実球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを算出し、かつ、画像中の各画素位置において、第１の半径より小さい第３の半径を有する第３の中実球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを算出し、第２の中実球を表す関数の１階偏微分ベクトルを固有ベクトルの方向に射影した値および第３の中実球を表す関数の１階偏微分ベクトルを固有ベクトルの方向に射影した値を用いて、第１の中実球を表す関数の各方向の２階偏微分の、第１の中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す補正を行う補正部を備えたものである。

本第２の発明による画像処理装置において、一方の位置の応答波形は、１つの正のピークと１つ負のピークが隣接する波形であり、第２の半径は、第１の中実球の中心から正のピークまでの長さと第１の中実球の中心から負のピークまでの長さのうち、長い方の長さと一致するものであり、第３の半径は、第１の中実球の中心から正のピークまでの長さと第１の中実球の中心から負のピークまでの長さのうち、短い方の長さと一致するものであることが好ましい。

上記、「第２の半径は、第１の中実球の中心から前記正のピークまでの長さと第１の中実球の中心から前記負のピークまでの長さのうち、長い方の長さと一致する」とは、第２の半径が第１の中実球の中心から正のピークまでの長さと第１の中実球の中心から負のピークまでの長さのうち、長い方の長さ（以下、第１の長さ）と厳密に一致する場合だけでなく、第２の中実球を表す関数の１階偏微分ベクトルを重み付けしたベクトルおよび第３の中実球を表す関数の１階偏微分ベクトルを重み付けしたベクトルを用いて、第１の中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す効果が得られる範囲であれば、第２の半径が第１の長さより大きい場合や小さい場合も含む。上記「一方の位置の応答波形を打ち消す効果」を良好に得るために、第２の半径と第１の長さができるだけ等しいことが好ましく、例えば、第２の半径と第１の長さとの差が２割以下であることが好ましく、１割以下であることがさらに好ましい。

同様に、上記、「第３の半径は、第１の中実球の中心から前記正のピークまでの長さと第１の中実球の中心から前記負のピークまでの長さのうち、短い方の長さと一致する」とは、第３の半径が第１の中実球の中心から正のピークまでの長さと第１の中実球の中心から負のピークまでの長さのうち、短い方の長さ（以下、第２の長さ）と厳密に一致する場合だけでなく、第２の中実球を表す関数の１階偏微分ベクトルを重み付けしたベクトルおよび第３の中実球を表す関数の１階偏微分ベクトルを重み付けしたベクトルを用いて、第１の中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す効果が得られる範囲であれば、第３の半径が第２の長さより大きい場合や小さい場合も含む。上記「一方の位置の応答波形を打ち消す効果」を良好に得るために、第３の半径と第２の長さができるだけ等しいことが好ましく、例えば、第３の半径と第２の長さとの差が２割以下であることが好ましく、１割以下であることがさらに好ましい。

本第１の発明による画像処理装置において、フィルタリング部は、複数のサイズの中実球を表す関数について、各中実球を表す関数の２階偏微分行列によるフィルタリングを施したヘッセ行列をそれぞれ算出するものであることが好ましい。

本第２の発明による画像処理装置において、フィルタリング部は、複数のサイズの第１の中実球を表す関数について、各第１の中実球を表す関数の２階偏微分行列によるフィルタリングを施したヘッセ行列をそれぞれ算出するものであることが好ましい。

本第１および第２の発明による画像処理装置において、画像が医用画像であり、構造物が血管であることが好ましい。

本第１の発明による画像処理装置において、フィルタリング部は、中実球を表す関数の２階偏微分行列を用いたフィルタリングを、フーリエ空間において実施するものであることが好ましい。また、本第２の発明による画像処理装置において、フィルタリング部は、第１の中実球を表す関数の２階偏微分行列を用いたフィルタリングを、フーリエ空間において実施するものであることが好ましい。

本第１および第２の発明による画像処理装置において、評価部は、構造物の局所的な点状構造、線状構造、面状構造の少なくとも１つを判別するものとすることが好ましい。

本第１発明によれば、画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、フィルタリング部が、画像中の各画素位置において、中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、取得した１階偏微分ベクトルを固有ベクトルの方向に射影した値を用いて、中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行う補正部を備えている。このため、中実球を表す関数が表す中実球の曲率より大きい曲率を有する構造物の輪郭部分と、中実球を表す関数の各方向の２階偏微分の応答波形の１つの位置のみが一致した場合に生じる誤判別を抑制でき、評価値の精度を向上することができる。このため、評価値に基づいて画像に含まれる構造物をより正確に判別することができる。

本第２発明によれば、画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、フィルタリング部が、第１の中実球の半径である第１の半径より大きい第２の半径を有する第２の中実球を表す関数の１階偏微分ベクトルを算出し、かつ、第１の半径より小さい第３の半径を有する第３の中実球を表す関数の１階偏微分ベクトルを算出し、第２の中実球を表す関数の１階偏微分ベクトルを固有ベクトルの方向に射影した値および第３の中実球を表す関数の１階偏微分ベクトルを固有ベクトルの方向に射影した値を用いて、中実球を表す関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す補正を行う補正部を備えている。このため、中実球を表す関数が表す中実球の曲率より大きい曲率を有する構造物の輪郭部分と、中実球を表す関数の各方向の２階偏微分の応答波形の１つの位置のみが一致した場合に生じる誤判別を抑制でき、評価値の精度を向上することができる。このため、評価値に基づいて画像に含まれる構造物をより正確に判別することができる。

本発明の第１の実施形態による画像処理装置の構成を示す概略ブロック図 多重解像度変換を説明するための図 線状構造の固有値を説明するための図 面状構造の固有値を説明するための図 本発明の第１の実施形態による補正処理の原理を説明するための図 本発明の第１の実施形態のフィルタリングに用いる中実球モデル関数のｘ方向の２階偏微分のレスポンスを説明するための図 本発明の第１の実施形態の補正処理に用いる中空球モデル関数のｘ方向の１階偏微分のレスポンスを説明するための図 本発明の第１の実施形態の中空球モデル関数の１階偏微分ベクトルにより補正された中実球のｘ方向の２階偏微分のレスポンスを説明するための図 本発明の第１の実施形態による補正処理を実施したレスポンスを説明する図（構造物のサイズと中実球モデル関数の表す中実球のサイズが一致する場合） 本発明の第１の実施形態による補正処理を実施したレスポンスを説明する図（構造物のサイズに対して中実球モデル関数の表す中実球のサイズが大きい場合） 本発明の第１の実施形態において行われる処理を示すフローチャート 本発明の第１の実施形態による画像処理を血管抽出に適用した例（擬似３次元画像） 本発明の第１の実施形態による画像処理を血管抽出に適用した例（断層像） 本発明の第２の実施形態による補正処理の原理を説明するための図 本発明の第２の実施形態の第２および第３の中実球モデル関数の１階偏微分ベクトルによる補正の原理を説明するための図

以下、図面を参照して本発明の実施形態について説明する。図１は本発明の実施形態による画像処理装置の構成を示す概略ブロック図である。なお、図１のような画像処理装置１の構成は、補助記憶装置（不図示）に読み込まれたプログラムをコンピュータ（例えばパーソナルコンピュータ等）上で実行することにより実現される。また、このプログラムは、ＣＤ−ＲＯＭ等の情報記憶媒体に記憶され、もしくはインターネット等のネットワークを介して配布され、コンピュータにインストールされることになる。

画像処理装置１は、例えばＸ線ＣＴ装置２により撮像された複数の２次元画像を用いて、３次元画像Ｍ０を生成し、この３次元画像Ｍ０に含まれる線状構造物と面状構造物とを自動的に分割するものであって、画像取得部１０、検出領域設定部２０、判別部３０、分割部４０、表示部５０および入力部６０を備える。

画像取得部１０は、例えばＸ線ＣＴ装置２により撮像された複数のＣＴ画像（２次元画像）を取得し、複数の２次元画像から３次元画像Ｍ０を生成する。なお、画像取得部１０は、ＣＴ画像のみならず、いわゆるＭＲＩ画像、ＲＩ画像、ＰＥＴ画像、Ｘ線画像等の２次元画像を取得するものであってもよい。

検出領域設定部２０は、まず３次元画像Ｍ０のボクセルサイズを等方化する。例えば、３次元画像Ｍ０のボクセルサイズが、３次元画像Ｍ０におけるＸ，Ｙ，Ｚ方向のそれぞれにおいて、０．３ｍｍ、０．３ｍｍ、０．６ｍｍであった場合、（Ｘ，Ｙ，Ｚ）＝（０．５，０．５，０．５）（ｍｍ）に等方化する。

検出領域設定部２０は、等方化後に３次元画像Ｍ０を多重解像度変換して図２に示すように、解像度が異なる複数の３次元多重解像度画像Ｍｓｉ（ｉ＝０〜ｎ）（ガウシアンピラミッド）を生成する。なお、ｉ＝０は３次元画像Ｍ０と同一解像度、ｉ＝ｎは最低解像度を表す。√２倍刻みで画像を縮小し、３次元多重解像度画像Ｍｓｉのボクセルサイズは、解像度が高い順に、（Ｘ，Ｙ，Ｚ）＝（０．５，０．５，０．５）、（０．７，０．７，０．７）、（１．０，１．０，１．０）…となる。

判別部３０は、フィルタリング部３２および補正部３３および評価部３４を備える。フィルタリング部３２は、ヘッセ行列（評価行列）を用いたヘッセ解析を行うために、３次元多重解像度画像Ｍｓｉのそれぞれに、後述の中実球を表す関数（中実球モデル関数）とガウシアンカーネルを用いた、非特許文献２による手法と同様のフィルタリングを行う。すなわち、解像度が異なる３次元多重解像度画像Ｍｓｉのそれぞれに対して、同一サイズのフィルタカーネルを畳み込む。なお、中実球モデル関数は中実球の半径Ｒとガウスカーネルのσで規定される。これらは、予備的な解析等の知見に基づいて応じた適切な値が設定される。なお、中実球の半径Ｒはガウスカーネルのσよりも少なくとも大きい値が設定される。

解像度が異なる３次元多重解像度画像Ｍｓｉのそれぞれに対して、同一サイズのフィルタカーネル（例えばR=2.0(voxel)、σ=0.5(voxel)とする）を畳み込むことにより、３次元画像Ｍ０に対して、実質的に異なるサイズのフィルタカーネルを適用することとなるため、異なるサイズの点状構造物、線状構造物（例えば、血管）および面状構造物（例えば皮質骨等の骨）を検出することができる。言い換えると、複数サイズの半径Ｒの中実球を表す関数を用いて、該各中実球を表す関数の２階偏微分行列を用いたフィルタリングを施した評価行列をそれぞれ算出している。

以下、非特許文献２による中実球モデル関数の２階偏微分を用いたヘッセ解析について説明する。ヘッセ解析に使用されるヘッセ行列は３次元画像に対しては下記の式（１）に示すように３×３の行列となる。

上記ヘッセ行列の各要素Ｉｘｘ、Ｉｘｙ、Ｉｙｙ、Ｉｙｚ、Ｉｚｚ、Ｉｚｘは、対象画像の画像データを、中空球モデル関数ｆ（ｒ）およびガウシアンカーネル関数ｇ（ｒ）の２階偏微分を用いてフィルタリング（コンボリューション演算）して算出する。

非特許文献２ではこのフィルタリング処理をフーリエ空間で行う。まず対象画像の画像データをフーリエ変換する。次いで、このフーリエ変換された対象画像（FT(image)）に、それぞれフーリエ変換された中実球モデル関数、ガウシアンカーネル関数、および２階偏微分を乗算する。そして、乗算結果を逆フーリエ変換することによってフィルタリング結果を取得する。なお、この取得されたフィルタリング結果は、実空間でコンボリューション演算をして得られたフィルタリング結果と同じものとなる。中実球モデル関数のフーリエ変換をF（ν）、ガウシアンカーネル関数のフーリエ変換をＧ（ν）としたとき、以上の関係は式（２）のように表される。

ここでｘ、ｙ、ｚは実空間における３次元座標であり、Ｒは中実球の半径である。また３次元フーリエ空間の変数（周波数）の極座標表現を式（４）に示したように、ν＝（ν_ｘ ^２＋ν_ｙ ^２＋ν_ｚ ^２）^１／２で表す。

式（２）に用いられている各関数について、フーリエ変換したガウシアンカーネル関数Ｇ（ν）を上記式（３）に示し、フーリエ変換した中実球モデル関数Ｆ（ν）を式（４）に示す。式（４）に示す中実球モデル関数Ｆ（ν）は、式（５）に示す中実球を表す関数ｆ（ｒ）をフーリエ変換することにより求められる。また、式（３）に示すガウシアンカーネル関数Ｇ（ν）は、非特許文献１および２と同様にフィルタリングの際に対象画像の画素値の微分範囲を規定するために用いられている。

また、式（２）において、（２πν_ｘ）^ｌ×（２πν_ｙ）^ｍ×（２πν_ｚ）^ｎの部分はフーリエ空間における微分処理に相当し、ｌ＋ｍ＋ｎ＝２（０＜ｌ、０＜ｍ、０＜ｎ）となるように、各微分方向に対応する係数ｌ、ｍ、ｎをそれぞれ式（４）に代入することにより、各要素Ｉｘｘ、Ｉｘｙ、Ｉｙｙ、Ｉｙｚ、Ｉｚｚ、Ｉｚｘに対応する値をそれぞれ算出できる。例えば、３次元多重解像度画像Ｍｓｉの処理の対象となる画素において、ヘッセ行列におけるＸ方向における２階偏微分である要素Ｉｘｘは、式（２）にＸ方向にｌ＝２を代入し、Ｙ方向およびＺ方向にそれぞれｍ＝ｎ＝０を代入することにより算出できる。また、３次元多重解像度画像Ｍｓｉの処理の対象となる画素において、ヘッセ行列における要素Ｉｘｙは、式（２）にＸ方向およびＹ方向のそれぞれにｌ＝ｍ＝１を代入し、Ｚ方向にｎ＝０を代入することにより算出できる。

このようにして算出したヘッセ行列を固有値分解して固有値を得たとき、線状構造は図３に示すように、３つのうち２つの固有値の値が大きく、１つが０に近い特徴を持つことが知られている。例えば、式（１）の固有値は、線状構造からなる対象組織に対して、式（６）に示す関係を有する。なお、固有値は｜λ_１｜＜｜λ_２｜＜｜λ_３｜とする。

また、面状構造は図４に示すように、３つのうち１つの固有値の値が大きく、２つが０に近い特徴を持つことが知られている。例えば、式（１）の固有値は、面状構造からなる対象組織に対して、式（７）のような関係を有する。

また、点状構造は、３つの固有値の全てが大きい特徴を持つことが知られている。例えば、式（１）の固有値は、点状構造からなる対象組織に対して、式（８）のような関係を有する。

したがって、固有値から線状構造らしさ、面状構造および点状構造らしさを判別することができ、判別結果を用いて、３次元画像Ｍ０において、線状構造物である血管領域、および面状構造物である骨領域を分割することができる。

本実施形態における補正部３３は、まず、フィルタリング部３２が算出したヘッセ行列を固有値分解し、３つの固有値λ_１，λ_２，λ_３（ただし、｜λ_１｜≦｜λ_２｜≦｜λ_３｜とする）を算出する。

ここで、図５Ａ、５Ｂ、５Ｃ、５Ｄを用いて本実施形態による誤判別の抑制の原理を説明し、その後、この原理に基づく具体的な補正部３３の補正処理を説明する。

図５Ａは、（Ａ）に中実球モデル関数（一点鎖線）と中実球モデル関数のｘ方向の２階偏微分（実線）の各ｘ位置のレスポンスを示し、（Ｂ）に中空球モデル関数（一点鎖線）と中空球モデル関数のｘ方向の１階偏微分（破線）の各ｘ位置のレスポンスを示し、（Ｃ）に中実球モデル関数と中空球モデル関数の１階偏微分を用いて補正された中実球モデル関数のｘ方向の２階偏微分（実線）の各ｘ位置のレスポンスを示す。

図５Ａ（Ａ）に示すように、中実球モデル関数の各方向の２階偏微分は中実球の表面に相当する離間した２つの位置に応答波形を示す。そして、血管の直径などの構造物を横断する線分と中実球モデル関数の表す中実球の直径が略一致した場合、各微分方向について画像中の構造物の輪郭のうち対向する２箇所の輪郭部分（構造物を横断する線分の両端）の両方が中実球モデル関数の２階偏微分の上記応答波形を示す２つの位置とそれぞれ一致して大きいレスポンス（期待したレスポンス）が得られるため、非特許文献２の手法によるヘッセ解析によれば血管などの構造物が判別できる。

一方、中実球より曲率の大きい構造物の輪郭部分において、画像中の構造物の輪郭部分に中実球モデル関数の各微分方向について２階偏微分による応答波形の１つの位置のみが一致した場合にも、この一つの位置の応答波形に対応する所定の大きさのレスポンスが得られる。本発明者は、この所定の大きさのレスポンスと期待したレスポンスとの区別がつかないことが中実球の直径と略同じ径の構造物が誤判別されてしまう原因となっていると推測した。そして、この中実球モデル関数の各方向の２階偏微分の１つの位置における応答波形を打ち消すことにより、誤判別が解消可能であることを見出した。そして、図５Ａ（Ａ）に示すように、中実球モデル関数の各方向２階偏微分は中実球の中心から対称な２つの位置に応答波形を示すという特徴を有するため、中実球モデル関数の２階偏微分の１つの応答波形と同じ位置にこの１つの応答波形と略同形状であるとともに正負の符号を逆にした応答波形を有する関数を利用することにより、中実球モデル関数の２階偏微分行列の１つの位置における応答波形を打ち消すことに注目した。

本実施形態では、図５Ａ（Ｂ）に示すように、中実球と略同じサイズの中空球モデル関数ｆ_１（ｒ）のｘ方向の１階偏微分のレスポンスが、中実球モデル関数のｘ方向の２階偏微分の一方の応答波形と同じ位置に、中実球モデル関数の１つの応答波形と略同形状かつ同符号の応答波形を有するという特徴を有し、さらに、中実球モデル関数のｘ方向の２階偏微分の他方の応答波形と同じ位置に、中実球モデル関数のｘ方向の２階偏微分の他方の応答波形と略同形状であるとともに正負の符号を逆にした応答波形を有するという特徴を有することを利用して、中実球と略同じサイズの中空球モデル関数ｆ_１（ｒ）ｘ方向の１階偏微分値を、中実球モデル関数のｘ方向の２階偏微分の１つの位置における応答波形を打ち消すために用いる。

つまり、中実球モデル関数のｘ方向の２階偏微分を用いてフィルタリングを行ったレスポンスと、中実球モデル関数の表す中実球と略同じ半径の中空球を表す中空球モデル関数の１階偏微分を用いてフィルタリングを行ったレスポンスとを加算すると、ｘ方向の負側の応答波形は、同じ位置において逆の符号を有するため打ち消し合い、ｘ方向の正側の応答波形は、同じ位置において同じ符号を有するため強め合うため、図５Ａ（Ｃ）に示すようにｘ方向において正側の１つの位置にのみ応答波形を有するものとなる。

図５Ｂ（Ａ）は、中実球モデル関数ｆ（ｒ）を示し、（Ｂ）は中実球モデル関数ｆ（ｒ）のｘ方向の２階偏微分のレスポンスを示す。図５Ｃ（Ａ）は、中空球モデル関数ｆ_１（ｒ）を示し、（Ｂ）は中空球モデル関数ｆ_１（ｒ）のｘ方向の１階偏微分のレスポンスを示す。図５Ｄは、図５Ｂの（Ｂ）に示す中実球モデル関数のｘ方向の２階偏微分のレスポンスと図５Ｃに示す中空球モデル関数ｆ_１（ｒ）の１階偏微分のレスポンスとを加算した結果を示す。なお、図５Ｂ〜Ｄにおいて、各上述のレスポンスを等間隔にｚ座標値を異ならせた複数のｘｙ平面図により示しており、これらの複数のｘｙ平面図をｚ座標の大きいものから垂直方向に上から順に並べて示している。

図５Ｂ〜Ｄの各レスポンスを表すｘｙ平面図において、明るい（白い）ほどレスポンスが正の方向に大きくなり、暗い（黒い）ほどレスポンスが負の方向に大きくなる。図５Ｂ（Ｂ）に示す中実球の中心から負側に位置する正のピークと負のピークの隣接した応答波形と図５Ｃ（Ｂ）に示す中実球の中心から負側に位置する負のピークと正のピークの隣接した応答波形が相互に打ち消し合い、図５Ｄでは、中実球の中心から負側には応答波形が見られなくなっているのが分かる。同様に各方向について中実球の中心から負側の応答波形を打ち消すことができる。

本実施形態では、上記原理に基づいて、補正部３３が中空球モデル関数の１階偏微分ベクトルを用いて中実球モデル関数ｆ（ｒ）の各方向の２階偏微分における１つの位置の応答波形を打ち消すようにヘッセ行列の固有値を補正する。具体的な補正方法を以下に説明する。

本実施形態において補正部３３は、まず、式（９）を用いて補正に用いるための１階偏微分ベクトルを算出する。詳細には、式（９）に示すように、フーリエ変換した３次元多重解像度画像Ｍｓｉの処理の対象となる画素に、フーリエ変換したガウシアンカーネル関数Ｇ（ν）と下記の式（１１）に表すフーリエ変換した中空球モデル関数Ｆ_１（ν）の１階偏微分フィルタの畳み込みを行い、そのフィルタリング結果を逆フーリエ変換することにより、ヘッセ行列の固有値補正に用いるための１階偏微分ベクトルを算出する。なお、式（１１）に用いた関数Ｆ_１（ν）は、式（１０）で表されるデルタ関数δ（ｒ−Ｒ_４）により定義される中空球モデル関数ｆ_１（ｒ）を、フーリエ変換することにより得られる。

また、式（９）において、（２πν_ｘ）^ｌ×（２πν_ｙ）^ｍ×（２πν_ｚ）^ｎの部分はフーリエ空間における微分処理に相当し、ｌ＋ｍ＋ｎ＝１（０＜ｌ、０＜ｍ、０＜ｎ）となるように、各微分方向に対応する係数ｌ、ｍ、ｎをそれぞれ式（９）に代入することにより、1階偏微分ベクトルの各要素ρ_１、ρ_２、ρ_３に対応する値をそれぞれ算出できる。

なお、式（９）においては、フーリエ空間においてフィルタリングを行っているが、実空間においてフィルタリングを行ってもよい。なお、式（１０）において、ｘ、ｙ、ｚは３次元空間の座標、ｒはその極座標表現であり、Ｒ_４は中空球の半径とする。また、本実施形態において中空球モデル関数の表す中空球の半径Ｒ_４は、中実球モデル関数の表す中実球の半径Ｒと同じである。なお、中空球の半径Ｒ_４は、上記中実球モデル関数の２階偏微分の１つの位置における応答波形を打ち消す効果を有する範囲で、中実球の半径Ｒよりも中空球の半径Ｒ_４が大きくてもよく、小さくてもよいが、中実球モデル関数の２階偏微分の応答波形と、中空球モデル関数の１階偏微分の応答波形をより対応したものとするために、中実球の半径Ｒと厳密に一致することが好ましい。

上述したように算出したＸ方向、Ｙ方向およびＺ方向の１階偏微分ベクトル（ρ_１,ρ_２,ρ_３）は、固有値λ_１，λ_２，λ_３の固有ベクトルｅ_１，ｅ_２，ｅ_３の方向とずれているため、補正部３３は、固有ベクトルｅ_１，ｅ_２，ｅ_３の方向に対応した１階偏微分ベクトルρ_１ ^’,ρ_２ ^’,ρ_３ ^’を下記の式（１２）により算出する。なお、評価行列の固有ベクトルをｅ_１＝（ｘ_１,ｙ_１、ｚ_１），ｅ_２＝（ｘ_２,ｙ_２、ｚ_２），ｅ_３＝（ｘ_３,ｙ_３、ｚ_３）とし、中空球を表す関数の１階偏微分ベクトルを（ρ_１,ρ_２,ρ_３）とする。

そして、ヘッセ行列の固有値を式（１３）に示すように補正する。なお、評価行列の固有値をλ₁、λ_２、λ_３とし、αを所定の係数とする。ここでαは中実球モデル関数の２階偏微分の原点（中実球の中心）から等間隔な２つの位置に表れる応答波形の一方を最も打ち消すように予め設計された重みである。この重みは中実球の半径Ｒとカーネルサイズsに合わせて設計しておく。

式（１３）において、λ_１＋αρ’は各微分方向において中実球モデル関数の中心に対称な２つの応答波形のうち、一方（例えば図５Ａ（Ａ）の左側の応答波形）を打ち消す補正を行ったときのレスポンスを表す。λ_１−αρ’はもう一方（例えば図５Ａ（Ａ）の右側の応答波形）を打ち消す補正を行ったときのレスポンスを表す。式（１３）に示すように、固有値λ_１’は、２つのレスポンス｜λ_１＋αρ’｜、｜λ_１−αρ’｜の小さいほうの値となるように補正されている。この結果、２つのレスポンス｜λ_１＋αρ’｜、｜λ_１−αρ’｜のいずれかの値が小さい場合には、各固有値λ_１’の値が小さくなり、２つのレスポンス｜λ_１＋αρ’｜、｜λ_１−αρ’｜がいずれも大きい値となる場合にのみ各固有値λ_１’の値が大きくなる。λ_２’、λ_３’ついても同様である。

従来の方法では誤って構造物を検出していた場合、つまり、中実球モデル関数の中心に対称な２つの応答波形の一方にのみレスポンスが得られる場合には、中実球モデル関数のいずれか一方の応答波形を打ち消すと、２つのレスポンス｜λ_１＋αρ’｜、｜λ_１−αρ’｜の少なくとも一方の値が小さくなるため、ｍｉｎ（｜λ_１＋αρ’｜、｜λ_１−αρ’｜）が小さい値となる。一方、構造物を正しく検出している場合、すなわち、各微分方向について画像中の構造物の輪郭のうち対向する２箇所の輪郭部分（構造物を横断する線分の両端）の両方が中実球モデル関数の２階偏微分の応答波形を示す２つの位置とそれぞれ一致している場合には、中実球モデル関数の中心に対称な２つの応答波形のうち、いずれか一方の応答波形を打ち消しても、もう一方の応答波形に基づいて大きいレスポンスが得られるため、２つのレスポンス｜λ_１＋αρ’｜、｜λ_１−αρ’｜の最小値が大きくなる。従って、誤って構造物を検出していた場合、つまり、２つの応答波形の一方にのみレスポンスが得られる場合を区別することができるため、正確に構造物を判別することができる。なお、式（１３）はレスポンスの絶対値を用いて補正をする例を示しているが、絶対値を用いないで補正を行ってもよい。また、ｍｉｎ（｜λ_１＋αρ’｜、｜λ_１−αρ’｜）を用いずにλ_１＋αρ’とλ_１−αρ’を別々に用いて対象構造の判別を行ってもよい。

そして、評価部３４は、下記の式（１４）、（１５）において、固有値λ_１，λ_２，λ_３に代えて補正後の固有値λ_１’，λ_２’，λ_３’を使用し、これにより算出した値Ｒ_A，Ｒ_B，Ｒ_Cを用いて３次元多重解像度画像Ｍｓｉの各画素における線状構造らしさの評価値Ｌ０（Lineness）および面状構造らしさの評価値Ｐ０（Planeness）を算出する。なお、先述の通り、評価部３４は、補正されたヘッセ行列の固有値λ_１’，λ_２’，λ_３’と固有ベクトルｅ_１，ｅ_２，ｅ_３に基づいて、点状構造、線状構造、面状構造を判別することが可能であるが、必ずしも点状構造、線状構造、面状構造の全てを評価する必要はなく、要求される仕様に応じて、点状構造、線状構造、面状構造の一部のみを判別してもよい。本実施形態では、医用画像から線状構造および面状構造を抽出することを目的としているため、線状構造および面状構造についてのみ評価値を算出している。

なお、式（１４）、（１５）におけるａ〜ｈは定数である。また、Ｒ_A，Ｒ_B，Ｒ_Cは下記の式（１６）〜（１９）により算出する。また、Ｓ_2ndは２階偏微分値のパワーであり、下記の式（１９）により算出する。

図６Ａ、６Ｂを用いて、線状構造物を表すサンプル画像に、フィルタリング部３２によるフィルタリング処理の後に補正部３３による補正処理を施したレスポンスの例を以下に示す。図６Ａは、構造物のサイズ（構造物を横断する線分の長さ）と中実球モデル関数のサイズ（中実球の直径）がほぼ一致する場合のレスポンスを説明する図である。図６Ａ（Ａ）は、線状構造物の輝度値のｘｙ平面図を示す。図６Ａ（Ｂ）は、図６Ａ（Ａ）の画像データに中実球モデル関数のｘ方向の２階偏微分を用いたフィルタリング処理を行って得られたレスポンスをｘｙ平面図で示す。図６Ａ（Ｃ）は図６Ａ（Ａ）の画像データに中空球モデル関数のｘ方向の１階偏微分を用いたフィルタリング処理を行って得られたレスポンスをｘｙ平面図で示す。図６Ａ（Ｄ）は、図６Ａ（Ａ）の画像に中実球モデル関数のｘ方向の２階偏微分行列を用いたフィルタリング処理を行って得られたレスポンスを中空球のｘ方向の１階偏微分のフィルタリング処理を行ったレスポンスを用いて補正したものをｘｙ平面図で表す。

図６Ａ（Ｄ）に示されるように、線状構造物の中心付近に線状構造らしさのピークが表れており、本実施形態による画像処理によって線状構造物を判別するための良好なレスポンスが得られていることが分かる。

図６Ｂは、ｘ方向の構造物のサイズ（構造物を横断する線分の長さ）が中実球モデル関数の中実球の直径よりも大きい場合のレスポンスを説明する図である。図６Ｂ（Ａ）は、線状構造物の輝度値のｘｙ平面図を示す。図６Ｂ（Ｂ）は、図６Ｂ（Ａ）画像データに中実球モデル関数のｘ方向の２階偏微分を用いたフィルタリング処理を行って得られたレスポンスをｘｙ平面図で示す。図６Ｂ（Ｃ）は、中空球モデル関数のｘ方向の１階偏微分を用いたフィルタリング処理を行って得られたレスポンスをｘｙ平面図で示す。図６Ｂ（Ｄ）は、図６Ｂ（Ａ）の画像データに中実球モデル関数のｘ方向の２階偏微分行列を用いたフィルタリング処理を行って得られたレスポンスを中空球のｘ方向の１階偏微分のフィルタリング処理を行ったレスポンスを用いて補正したものをｘｙ平面図で表す。

非特許文献２の手法に対応する評価値を示す図６Ｂ（Ｂ）では中実球モデル関数の２階偏微分による誤判別の原因となる一定のレスポンスが２箇所に表れている（２つの負のピーク）のに対し、本実施形態の手法によるレスポンスを示す図６Ｂ（Ｄ）では、レスポンスが非常に弱くなり、誤判別の原因となるレスポンスが抑制されていることが分かる。

なお、上記図６Ａ（Ｄ）、図６Ｂ（Ｄ）は、式（１３）のλ_１に代えて、式（１）に示すヘッセ行列の要素（Ｉｘｘ）を適用した場合のレスポンスを示している。固有値を算出する処理は、ヘッセ行列における空間を線形変換しているだけでレスポンスの値に影響を与えるものではないため、上記式（１３）において、ヘッセ行列を固有値解析して得られた固有値λ_１に替えて、固有値解析をする前のヘッセ行列の要素を補正しても、補正された要素において誤検出の原因となるレスポンスが打ち消されている効果が観察できるからである。

なお、本実施形態においては、解像度が異なる３次元多重解像度画像Ｍｓｉにおいて、それぞれ線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０が算出される。算出された評価値Ｌ０，Ｐ０は、元の３次元画像Ｍ０の対応する画素位置の評価値となるが、すべての３次元多重解像度画像Ｍｓｉの対応する画素位置において評価値が算出されることとなる。

分割部４０は、判別部３０が算出した線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０に基づいて、３次元画像Ｍ０の血管領域および血管以外の領域を領域分割する。具体的には、血管領域を対象領域、血管領域以外の領域を背景領域に設定し、３次元画像Ｍ０内の全画素位置において所定画素サイズの判別領域を設定し、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領域分割方法を用いて、３次元画像Ｍ０を対象領域（血管領域）と背景領域とに分割する（ステップＳＴ６）。そして、表示部５０が分割された対象領域および背景領域をボリュームレンダリング表示し（ステップＳＴ７）、処理を終了する。

本実施形態の画像処理を実際の患者の胸部ＣＴ画像の血管判別に適用した例と従来の方法を同一のＣＴ画像の血管判別にそれぞれ適用した例を図８Ａ、８Ｂに示す。図８Ａ、８Ｂにおいて、左が従来手法による画像処理の適用例（比較例）であり、右が本実施形態の画像処理の適用例である。図８Ａは、判別された血管領域をボリュームレンダリング法により表した擬似３次元画像を示し、図８Ｂは、図８Ａの一部を拡大したアキシャル断層像で表したものである。

図８Ａに示すように従来の手法では、多数の誤判別による血管が表示されているところ、本実施形態の適用例においては大幅に誤判別が抑制されて正確に血管が抽出されていることが分かる。また、図８Ｂの比較例（左図）では、血管の内側に線状構造物らしいと誤って判別された部分が薄い灰色で示されているが、本実施形態の適用例（右図）では、血管の内側に線状構造物らしいと判断された部分（薄い灰色の部分）がないことが分かる。

このように、本実施形態によれば、中実球モデル関数の２階偏微分の１つの応答波形と同じ位置にこの１つの応答波形と略同形状であるとともに正負の符号を逆にした応答波形を有する関数の１階偏微分ベクトルを算出し、この１階偏微分ベクトルをヘッセ行列の固有ベクトルの方向に射影した値を用いて、中実球モデル関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち、１つの位置に表れる応答波形を打ち消すように補正したため、構造物の輪郭と中実球モデル関数のいずれかの方向の２階偏微分の応答波形の表れる１つの位置のみが一致した場合に生じる誤判別を抑制でき、評価値の精度を向上することができる。このため、評価値に基づいて画像に含まれる構造物をより正確に判別することができる。一般にフーリエ空間における画像のフィルタリング処理において、フィルタ係数の一部の好ましくない値のみを変更する補正を行うことは容易でないが、上記本実施形態のように中実球モデル関数の２階偏微分のレスポンス特性に着目して中空球モデル関数の１階偏微分のレスポンス特性を利用することにより、良好に非特許文献１の手法による誤判別の問題を解消することができる。

また、第１の実施形態においては中空球モデル関数の「１階」偏微分ベクトルを補正に用いているため、一旦ｘ、ｙ、ｚ方向の偏微分値を計算すると、任意の方向の偏微分値を計算することができる。すなわち、中空球モデル関数の「１階」偏微分ベクトルを２階偏微分の固有ベクトルの方向に射影して用いることにより、任意の方向の２階偏微分の補正が可能である。結果として、補正のための計算処理の増加が比較的小さく、様々な処理能力の計算機を本実施形態の画像処理装置として使用可能である。

また、中空球モデル関数の各方向の１階偏微分の１つの位置の応答波形が、中実球モデル関数の各方向の２階偏微分の１つの位置の応答波形と、正負の符号が逆であるとともに非常に近似した形状であるため、誤判別の原因となる中実球モデル関数の２階偏微分の１つの位置の応答波形を好適に打ち消すことができる。

また、上記第１の実施形態においては、複数サイズの半径Ｒの中実球を表す関数について、ヘッセ行列をそれぞれ算出し、この複数のヘッセ行列に基づいて、それぞれ誤判別の原因となるレスポンスを抑制する補正を行って、各画素位置における線状構造らしさおよび面状構造らしさの評価値を算出している。このため、各中実球サイズに対応するサイズの構造物を好適に評価することができる。また、この評価値に基づいて、上記の第１の実施形態のように対象領域を分割した場合には、好適に複数サイズの構造物を対象領域として分割することができる。

なお、異なるサイズの構造物に対してヘッセ解析を行う手法として、本実施形態のように複数の解像度の画像に対して、一定のカーネルサイズσおよび一定の中実球モデル関数の半径Ｒを適用してヘッセ解析を行う手法を適用してもよく、一定の解像度の画像に対して、カーネルサイズσおよび中実球モデル関数の半径Ｒを異ならせてヘッセ解析を行う手法を適用してもよい。例えば、構造物の長軸や短軸の長さなど、検出対象構造物を横断する線分の長さを試験的に得られたデータなどから予め取得しておき、検出対象構造物を横断する線分の長さと半径Ｒの大きさが一致するように、３次元多重解像度画像Ｍｓｉの解像度、あるいは、カーネルサイズσおよび中実球の半径Ｒを設定することが考えられる。

また、上記第１の実施形態においては、線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０の双方を算出しているが、線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０のいずれか一方のみを算出するようにしてもよい。

本発明の第２の実施形態を以下に説明する。第２の実施形態においては、補正部３３による補正処理が第１の実施形態と異なるが、それ以外は第１の実施形態と同じであり、各機能ブロックの機能も共通している。以下、第１の実施形態と異なる点を中心に説明し、第１の実施形態と同じ点については説明を省略する。

先述の中実球モデル関数の２階偏微分による２つの位置における応答波形の１つの位置における応答波形を打ち消すために、第２の実施形態では、中実球モデル関数の２階偏微分のｘ方向の負側に位置する応答波形がｘ方向の負側から１つの正のピークとこれに隣接する１つの負のピークを有するものであり、ｘ方向の正側に位置する応答波形がｘ方向の負側から１つの負のピークとこれに隣接する１つの正のピークを有するものであるという特徴に着目し、この一方の応答波形と同じ位置に位置し、かつ、この一方の応答波形と正負の符号が逆で略同形状の応答波形を、この中実球よりやや小さいサイズの中実球モデル関数ｆ_２（ｒ）（第２の中実球モデル関数）の各方向の１階偏微分と、この中実球よりやや大きいサイズの中実球モデル関数ｆ_３（ｒ）（第３の中実球モデル関数）の各方向の１階偏微分を組み合わせて作成可能であることを見出した。なお、ここでは、区別のために、フィルタリング部３２がフィルタリング処理に用いる中実球モデル関数を第１の中実球モデル関数と称する。

図９Ａ（Ａ）は、第２の半径Ｒ_２を有する第２の中実球モデル関数ｆ_２（ｒ）のｘ方向の１階偏微分のレスポンスを示し、図９Ａ（Ｂ）は、第３の半径Ｒ_３を有する第３の中実球モデル関数ｆ_３（ｒ）のｘ方向の１階偏微分のレスポンスを示し、図９Ａ（Ｃ）は、正負の符号を反転させた第２の中実球モデル関数ｆ_２（ｒ）のｘ方向の１階偏微分に第３の中実球モデル関数ｆ_３（ｒ）のｘ方向の１階偏微分を加算したレスポンスを示す。

図９Ａ（Ａ）、（Ｂ）に示すように、第２および第３の中実球モデル関数ｆ_２（ｒ）、ｆ_３（ｒ）の１階偏微分のレスポンスは、中実球の中心から負側の位置に１つの正のピークを有する応答波形を有し、中実球の中心から正側の位置に１つの負のピークを有する応答波形を有し、各応答波形は中心から対称に表れるという特徴を有する。

そして、図９Ａ（Ａ）に示す中実球モデル関数ｆ_２（ｒ）のｘ方向の１階偏微分の正負の符号を反転させたものと、図９Ａ（Ｂ）に示す中実球モデル関数ｆ_３（ｒ）のｘ方向の１階偏微分を組み合わせて用いることにより、図９Ａ（Ｃ）に示すようなｘ方向のレスポンスを実現できる。図９Ａ（Ｃ）は、中実球の中心から対称な２つの位置に応答波形を有するものであり、ｘ方向の負側に位置する応答波形はｘ方向の負側から１つの負のピークとこれに隣接する１つの正のピークを有するものであり、ｘ方向の正側に位置する応答波形はｘ方向の負側から１つの負のピークとこれに隣接する１つの正のピークを有するという特徴を有するものである。

つまり、中実球モデル関数の各方向の２階偏微分行列によるフィルタリングに、第２および第３の中実球モデル関数ｆ_２（ｒ）、ｆ_３（ｒ）の各方向の１階偏微分を図９Ａ（Ｃ）のようなレスポンスを組み合わせると、各方向の正側の応答波形は、同じ位置において逆の符号を有するため打ち消し合い、ｘ方向の負側の応答波形は、同じ位置において同じ符号を有するため強め合うため、図５Ａ（Ｃ）に示すようなｘ方向の正側の１つの位置にのみ応答波形を有するものとなる。

なお、本実施形態では、第１の中実球モデル関数の２階偏微分の１つの位置の応答波形の打ち消しを行うために、第２の半径Ｒ_２は、第１の中実球モデル関数の２階偏微分における中実球の中心から正方向に位置する応答波形の正のピークまでの距離と同じであり、第３の半径Ｒ_３は、第１の中実球モデル関数の２階偏微分における中実球の中心から正方向に位置する応答波形の負のピークまでの距離と同じである。この第２の半径Ｒ_２（第３の半径Ｒ_３）は、上記中実球モデル関数の各方向の２階偏微分の１つの位置における応答波形を打ち消す効果を有する範囲で、第２の半径Ｒ_２（第３の半径Ｒ_３）より「中心から応答波形の正（負）のピークまでの距離」が大きくてもよく、小さくてもよいが、中実モデル関数の２階偏微分の応答波形と、第２および第３の中実球モデル関数ｆ_２（ｒ）、ｆ_３（ｒ）の各方向の１階偏微分を組み合わせた応答波形をより対応したものとするために、本実施形態のように第２の半径Ｒ_２（第３の半径Ｒ_３）が「中心から応答波形の正（負）のピークまでの距離と厳密に一致することが好ましい。

ｘ方向における第２および第３の中実球モデル関数の１階偏微分のレスポンスを図９Ｂを用いて説明する。図９Ｂ（Ａ）は、第２の中実球モデル関数ｆ２（ｒ）のｘ方向の１階偏微分のレスポンスを示し、図９Ｂ（Ｂ）は第３の中実球モデル関数ｆ３（ｒ）のｘ方向の１階偏微分のレスポンスを示し、図９Ｂ（Ｃ）は、正負の符号を反転させた中実球モデル関数ｆ２（ｒ）のｘ方向の１階偏微分に中実球モデル関数ｆ３（ｒ）のｘ方向の１階偏微分を加算した場合のレスポンスを示す。なお、図９Ｂにおいて、各上述のレスポンスを等間隔にｚ座標値を異ならせた複数のｘｙ平面図により示しており、これらの複数のｘｙ平面図をｚ座標の大きいものから垂直方向に並べて示している。図９Ｂにおける各ｘｙ平面図において、明るい（白い）ほどレスポンスが正の方向に大きくなり、暗い（黒い）ほどレスポンスが負の方向に大きくなる。

図９Ｂ（Ｃ）に示す負側の応答波形は、図９Ｂ（Ａ）に示す第２の中実球モデル関数のｘ方向の１階偏微分の中心から負側に位置する負のピークと、図９Ｂ（Ｂ）に示す第３の中実球モデル関数のｘ方向の１階偏微分の中心から負側に位置する負のピークがずれて重なり合い、負の方向から負のピークと正のピークが隣接する波形となっていることが分かる。また、図９Ｂ（Ｃ）に示す正側の応答波形は、図９Ｂ（Ａ）に示す第２の中実球モデル関数の１階偏微分の中心から正側に位置する正のピークと、図９Ｂ（Ｂ）に示す第３の中実球モデル関数の１階偏微分の中心から正側に位置する正のピークがずれて重なり合い、負の方向から負のピークと正のピークが隣接する波形となっていることが分かる。第２および第３の中実球モデル関数の１階偏微分を図９Ｂ（Ｃ）のレスポンスを示すように組み合わせて、図５Ｃ（Ｂ）に示す第１の中実球モデル関数のｘ方向の２階偏微分を補正すると、図５Ｄに示すように中実球の中心から負側の応答波形を打ち消すことができる。同様に各方向について中実球の中心から負側の応答波形を打ち消すことができる。

本第２の実施形態では、上記原理に基づいて、補正部３３は、第２および第３の中実球モデル関数の各１階偏微分を用いて第１の中実球モデル関数ｆ（ｒ）の１つの位置の応答波形のみを打ち消すようにヘッセ行列を補正する。具体的な補正方法を以下に説明する。

まず、補正部３３は、第２の半径Ｒ_２の中実球モデル関数について、式（２）において半径Ｒに替えて第２の半径Ｒ_２を用いることにより第２の中実球モデル関数の１階偏微分ベクトルを算出する。なお、式（２）では、（２πν_ｘ）^ｌ×（２πν_ｙ）^ｍ×（２πν _ｚ ）^ｎの部分はフーリエ空間における微分処理に相当し、ｌ＋ｍ＋ｎ＝１（０＜ｌ、０＜ｍ、０＜ｎ）となるように、各微分方向に対応する係数ｌ、ｍ、ｎをそれぞれ式（９）に代入することにより、第２の中実球モデル関数ｆ_２（ｒ）について算出した各１階偏微分ベクトル（ｓ_ｘ、ｓ_ｙ、ｓ_ｚ）をそれぞれ算出できる。また、同様に、補正部３３は、第３の中実球モデル関数ｆ_３（ｒ）についても、式（２）において半径Ｒに替えて第３の半径Ｒ３を用いることにより１階偏微分ベクトル（ｔ_ｘ、ｔ_ｙ、ｔ_ｚ）を算出する。そして、第２の中実球モデル関数ｆ_２（ｒ）について算出した各１階偏微分ベクトル（ｓ_ｘ、ｓ_ｙ、ｓ_ｚ）、および第３の中実球モデル関数ｆ_３（ｒ）について算出した各１階偏微分ベクトル（ｔ_ｘ、ｔ_ｙ、ｔ_ｚ）を重み付け加算して、（ρ_１、ρ_２、ρ_３）＝（ｓ_ｘ−βｔ_ｘ、ｓ_ｙ−βｔ_ｙ、ｓ_ｚ−βｔ_ｚ）を算出する。

なお、式（２）においては、フーリエ空間においてフィルタリングを行っているが、実空間においてフィルタリングを行ってもよい。

上述したように算出したＸ方向、Ｙ方向およびＺ方向の１階偏微分ベクトル（ρ_１,ρ_２,ρ_３）は、固有値λ_１，λ_２，λ_３の固有ベクトルｅ_１，ｅ_２，ｅ_３の方向とずれているため、補正部３３は、固有ベクトルｅ_１，ｅ_２，ｅ_３の方向に対応した１階偏微分ベクトルρ_１ ^’,ρ_２ ^’,ρ_３ ^’を下記の式（１２）により算出する。

そして、第１の実施形態同様に、式（１３）に基づいてヘッセ行列の固有値を補正する。なお、評価行列の固有値をλ₁、λ_２、λ_３とし、α、βを所定の係数とする。ここでα、βは中実球モデル関数の２階偏微分の原点（中実球の中心）から等間隔な２つの位置に表れる応答波形の一方を最も打ち消すように予め設計された重みである。この重みは中実球の半径Ｒとカーネルサイズsに合わせて設計しておく。

そして、評価部３４は、第１の実施形態と同様に、式（１４）、（１５）において、固有値λ_１，λ_２，λ_３に代えて固有値λ_１’，λ_２’，λ_３’を使用し、これにより算出した値Ｒ_A，Ｒ_B，Ｒ_Cを用いて３次元多重解像度画像Ｍｓｉの各画素における線状構造らしさの評価値Ｌ０（Lineness）および面状構造らしさの評価値Ｐ０（Planeness）を算出する。

第２の実施形態においても、第１の実施形態同様に、中実球モデル関数の２階偏微分の１つの応答波形と同じ位置にこの１つの応答波形と略同形状であるとともに正負の符号を逆にした応答波形を有する関数（２種類のサイズの中実球モデル関数の１階偏微分を組み合わせた関数）の１階偏微分ベクトルを算出し、この１階偏微分ベクトルをヘッセ行列の固有ベクトルの方向に射影した値を用いて、中実球モデル関数の各方向の２階偏微分の、中実球の中心から対称に離間した２つの位置に表れる応答波形のうち、１つの位置に表れる応答波形を打ち消すように補正している。このため、第２の実施形態では、２種類のサイズの中実球モデル関数のそれぞれの１階偏微分を補正に用いる分、１種類のサイズの中空球モデル関数の１階偏微分を補正に用いる第１の実施形態の手法よりも補正部３３による補正処理の計算量は若干増加するものの、第１の実施形態同様に構造物の輪郭と中実球モデル関数のいずれかの方向の２階偏微分の応答波形の表れる１つの位置のみが一致した場合に生じる誤判別を抑制でき、評価値の精度を向上することができる。

また、第２の実施形態において、第２の半径Ｒ_２は、中実球の中心から正のピークまでの長さと中実球の中心から負のピークまでの長さのうち、長い方の長さと一致するものであり、第３の半径Ｒ_３は、中実球の中心から正のピークまでの長さと中実球の中心から負のピークまでの長さのうち、短い方の長さと一致するものであるため、応答波形を打ち消す効果が、第１の中実球の正のピークと負のピークが隣接した一方の応答波形により対応したものとなり、好適に誤判別の原因となる応答波形の抑制を実現できる。

また、上記第２の実施形態においても、複数サイズの半径Ｒの第１の中実球を表す関数について、ヘッセ行列をそれぞれ算出し、この複数のヘッセ行列に基づいて、それぞれ誤判別の原因となるレスポンスを抑制する補正を行って、各画素位置における線状構造らしさおよび面状構造らしさの評価値を算出することが好ましい。この場合には、各カーネルサイズに対応するサイズの構造物を好適に評価することができる。またこのために、複数の解像度の画像に対して一定の半径Ｒの第１の中実球モデル関数を用いて第２の実施形態の画像処理を行ってもよく、１種類の解像度の画像に対して、第１の半径Ｒを異ならせて、複数サイズの第１の半径Ｒの第１の中実球モデル関数を用いて第２の実施形態の画像処理を行ってもよい。

上記各実施形態による画像処理方法は、以下のような方向判別のための画像処理装置にも適用可能である。例えば、方向判別のための画像処理装置は、第１の実施形態における分割部４０と評価部３４と、表示部５０を備えておらず、代わりに方向判別部を備えたものとできる。このような画像処理装置の一例として、画像取得部１０による画像取得処理と、検出領域設定部２０による検出領域設定処理と、フィルタリング部３２によるフィルタリング処理については、第１の実施形態と同じであり、各機能ブロックの機能も共通したものとして、方向判別部を以下のように構成できる。

方向判別部は、フィルタリング部３２が算出した評価行列を固有値解析して固有値および固有ベクトルｅ_１，ｅ_２，ｅ_３を算出し、算出された固有ベクトルの向きに基づいて線状構造物の主軸方向または面状構造の法線方向を判別する。詳細には、フィルタリング部３２が得たヘッセ行列の固有値λ_１，λ_２，λ_３および固有ベクトルｅ_１，ｅ_２，ｅ_３を取得し、式（６）に示す関係を満たしていれば、固有値λ_１に対応する固有ベクトルｅ_１の指す方向を線状構造の軸方向と判別し、式（７）に示す関係を満たしていれば、固有値λ_３に対応する固有ベクトルｅ_３の指す方向を面状構造の法線方向と判別する。なお、方向判別部は、ヘッセ行列の固有値と固有ベクトルに基づいて、式（６）に基づく線状構造の軸方向、または、式（７）に基づく面状構造の法線方向のいずれか一方のみを判別してもよい。

上記のような方向を判別するための画像処理装置により得られた線状構造の軸方向や面状構造の法線方向は、例えば血管や大腸など管状構造のＣＰＲ(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アルゴリズムを実行することによって、線状構造物および面状構造物を分割することが可能である。

上記の各実施形態はあくまでも例示であり、上記のすべての説明が本発明の技術的範囲を限定的に解釈するために利用されるべきものではない。

この他、上記の実施形態におけるシステム構成、ハードウェア構成、処理フロー、モジュール構成、ユーザインターフェースや具体的処理内容等に対して、本発明の趣旨から逸脱しない範囲で様々な改変を行ったものも、本発明の技術的範囲に含まれる。

１ 画像処理装置
２ Ｘ線ＣＴ装置
１０ 画像取得部
２０ 検出領域設定部
３０ 判別部
３２ フィルタリング部
３３ 補正部
３４ 評価部
４０ 分割部
５０ 表示部
６０ 入力部

Claims (13)

1. 画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、
   前記算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、
   前記フィルタリング部が、前記画像中の前記各画素位置において、前記中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、前記取得した１階偏微分ベクトルを前記固有ベクトルの方向に射影した値を用いて、前記中実球を表す関数の各方向の２階偏微分の、前記中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行う補正部を備えたものであることを特徴とする画像処理装置。
2. 前記ヘッセ行列の固有値をλ₁、λ_２、λ_３とし、固有ベクトルｅ_１＝（ｘ_１,ｙ_１ ,ｚ_１），ｅ_２＝（ｘ_２,ｙ_２ ,ｚ_２），ｅ_３＝（ｘ_３,ｙ_３ ,ｚ_３）とし、前記中空球を表す関数の前記１階偏微分ベクトルを（ρ_１,ρ_２,ρ_３）とすると、前記補正部は、下記式（１２）で算出されるρ_１ ^’,ρ_２ ^’,ρ_３ ^’および所定の係数αを用いて、下記式（１３）に示すように、前記固有値を補正することにより、前記一方の応答波形を打ち消す補正を行うものであることを特徴とする請求項１記載の画像処理装置。
   
3. 前記中空球を表す関数は、下記式（１０）により表されるものであることを特徴とする請求項１または２記載の画像処理装置。
   ここで、ｘ、ｙ、ｚは３次元空間の座標、ｒはその極座標表現であり、Ｒ_４は前記中空球の半径とする。
4. 前記フィルタリング部は、複数のサイズの前記中実球を表す関数について、該各中実球を表す関数の２階偏微分行列によるフィルタリングを施した前記ヘッセ行列をそれぞれ算出するものであることを特徴とする請求項１から３のいずれか１項記載の画像処理装置。
5. 前記フィルタリング部は、前記中実球を表す関数の２階偏微分行列を用いたフィルタリングを、フーリエ空間において実施するものであることを特徴とする請求項１から４のいずれか１項記載の画像処理装置。
6. 画像中の各画素位置において、第１の中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング部と、
   前記算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価部とを備え、
   前記フィルタリング部が、前記画像中の各画素位置において、前記第１の中実球の半径である第１の半径より大きい第２の半径を有する第２の中実球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを算出し、かつ、前記画像中の各画素位置において、前記第１の半径より小さい第３の半径を有する第３の中実球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを算出し、前記第２の中実球を表す関数の１階偏微分ベクトルを前記固有ベクトルの方向に射影した値および第３の中実球を表す関数の１階偏微分ベクトルを前記固有ベクトルの方向に射影した値を用いて、前記第１の中実球を表す関数の各方向の２階偏微分の、前記第１の中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の位置の応答波形を打ち消す補正を行う補正部を備えたものであることを特徴とする画像処理装置。
7. 前記一方の位置の応答波形は、１つの正のピークと１つ負のピークが隣接する波形であり、
   前記第２の半径は、前記第１の中実球の中心から前記正のピークまでの長さと前記第１の中実球の中心から前記負のピークまでの長さのうち、長い方の長さと一致するものであり、
   前記第３の半径は、前記第１の中実球の中心から前記正のピークまでの長さと前記第１の中実球の中心から前記負のピークまでの長さのうち、短い方の長さと一致するものであることを特徴とする請求項６に記載の画像処理装置。
8. 前記フィルタリング部は、複数のサイズの前記第１の中実球を表す関数について、該各第１の中実球を表す関数の２階偏微分行列によるフィルタリングを施した前記ヘッセ行列をそれぞれ算出するものであることを特徴とする請求項６または７のいずれか１項記載の画像処理装置。
9. 前記フィルタリング部は、前記第１の中実球を表す関数の２階偏微分行列を用いたフィルタリングを、フーリエ空間において実施するものであることを特徴とする請求項６から８のいずれか１項記載の画像処理装置。
10. 前記評価部は、前記構造物の局所的な点状構造、線状構造、面状構造の少なくとも１つを判別するものであることを特徴とする請求項１から９のいずれか1項記載の画像処理装置。
11. 前記画像が医用画像であり、前記構造物が血管であることを特徴とする請求項１から１０のいずれか１項記載の画像処理装置。
12. 画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング工程と、
    前記算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価工程とを備えた画像処理方法であって、
    前記フィルタリング工程が、前記画像中の前記各画素位置において、前記中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、前記取得した１階偏微分ベクトルを前記固有ベクトルの方向に射影した値を用いて、前記中実球を表す関数の各方向の２階偏微分の、前記中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行うものであることを特徴とする画像処理方法。
13. コンピュータに、
    画像中の各画素位置において、中実球を表す関数の２階偏微分を用いたフィルタリングを施してヘッセ行列を算出するフィルタリング工程と、
    前記算出されたヘッセ行列を固有値解析して得られた固有値および固有ベクトルを用いて、画像中に含まれる構造物を判別する評価工程とを実行させる画像処理プログラムであって、
    前記フィルタリング工程が、前記画像中の前記各画素位置において、前記中実球の半径と同じ半径を有する中空球を表す関数の１階偏微分を用いたフィルタリングを施して１階偏微分ベクトルを取得し、前記取得した１階偏微分ベクトルを前記固有ベクトルの方向に射影した値を用いて、前記中実球を表す関数の各方向の２階偏微分の、前記中実球の中心から対称に離間した２つの位置に表れる応答波形のうち一方の応答波形を打ち消す補正を行うものであることを特徴とする画像処理プログラム。