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

画像処理装置および方法並びにプログラム Download PDF

Info

Publication number
JP5748693B2
JP5748693B2 JP2012055599A JP2012055599A JP5748693B2 JP 5748693 B2 JP5748693 B2 JP 5748693B2 JP 2012055599 A JP2012055599 A JP 2012055599A JP 2012055599 A JP2012055599 A JP 2012055599A JP 5748693 B2 JP5748693 B2 JP 5748693B2
Authority
JP
Japan
Prior art keywords
evaluation
image
filtering
matrix
hollow
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2012055599A
Other languages
English (en)
Other versions
JP2013188289A (ja
Inventor
嘉郎 北村
嘉郎 北村
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fujifilm Corp
Original Assignee
Fujifilm Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Fujifilm Corp filed Critical Fujifilm Corp
Priority to JP2012055599A priority Critical patent/JP5748693B2/ja
Priority to PCT/JP2013/001528 priority patent/WO2013136750A1/ja
Publication of JP2013188289A publication Critical patent/JP2013188289A/ja
Priority to US14/481,741 priority patent/US9307948B2/en
Application granted granted Critical
Publication of JP5748693B2 publication Critical patent/JP5748693B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5211—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
    • A61B6/5217—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data extracting a diagnostic or physiological parameter from medical diagnostic data
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/50—Clinical applications
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/50—Clinical applications
    • A61B6/504—Clinical applications involving diagnosis of blood vessels, e.g. by angiography
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/13—Edge detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/162—Segmentation; Edge detection involving graph-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/194—Segmentation; Edge detection involving foreground-background segmentation
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • G16H50/30—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for calculating health indices; for individual health risk assessment
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/02—Devices for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computerised tomographs
    • A61B6/032—Transmission computed tomography [CT]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/50—Clinical applications
    • A61B6/505—Clinical applications involving diagnosis of bone
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20024—Filtering details
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30061—Lung

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)）に、それぞれフーリエ変換された中空球モデル関数、ガウシアンカーネル関数、および２階偏微分を積算する。そして、積算結果を逆フーリエ変換することによってフィルタリング結果を取得する。なお、この取得されたフィルタリング結果は、実空間でコンボリューション演算をして得られたフィルタリング結果と同じものとなる。中空球モデル関数のフーリエ変換をＦ（ν）、ガウシアンカーネル関数のフーリエ変換をＧ（ν）としたとき、以上の関係は式（２）のように表される。
ここで、本実施形態における中空球モデル関数ｆ（ｒ）は、式（３）に示すようにデルタ関数δ（ｒ−Ｒ）により定義される。なお、フーリエ空間においてフィルタリングを行うため式（３）をフーリエ変換した式（４）が式（２）において用いられる。
ただし、式（３）のｘ、ｙ、ｚは実空間における３次元座標であり、Ｒは中空球の半径である。また３次元フーリエ空間の変数（周波数）の極座標表現をν＝（νｘ ２＋νｙ ２＋νｚ ２）１／２で表す。
また、式（２）におけるフーリエ変換されたガウシアンカーネル関数Ｇ(ν)は次式で表される。
なお、このガウシアンカーネル関数Ｇ（ν）は、非特許文献１および２と同様にフィルタリングの際に対象画像の画素値の微分範囲に幅を持たせるために用いる。
また、式（２）において、（２πνｘ）ｌ×（２πνｙ）ｍ×（２πνｚ）ｎの部分はフーリエ空間における微分処理に相当し、ｌ＋ｍ＋ｎ＝２（０＜ｌ、０＜ｍ、０＜ｎ）となるように、各微分方向に対応する係数ｌ、ｍ、ｎをそれぞれ式（２）に代入することにより、各要素Ｉｘｘ、Ｉｘｙ、Ｉｙｙ、Ｉｙｚ、Ｉｚｚ、Ｉｚｘに対応する値をそれぞれ算出できる。例えば、３次元多重解像度画像Ｍｓｉの処理の対象となる画素において、ヘッセ行列におけるＸ方向における二階偏微分値である要素Ｉｘｘは、式（２）にＸ方向にｌ＝２を代入し、Ｙ方向およびＺ方向にそれぞれｍ＝ｎ＝０を代入することにより算出できる。また、３次元多重解像度画像Ｍｓｉの処理の対象となる画素において、ヘッセ行列における要素Ｉｘｙは、式（２）にＸ方向およびＹ方向のそれぞれにｌ＝ｍ＝１を代入し、Ｚ方向にｎ＝０を代入することにより算出できる。
このようにして算出したヘッセ行列を固有値分解して固有値を得たとき、線状構造は図４に示すように、３つのうち２つの固有値の値が大きく、１つが０に近い特徴を持つことが知られている。例えば、式（１）の固有値は、線状構造からなる対象組織に対して、式（６）に示す関係を有する。なお、固有値は｜λ１｜＜｜λ２｜＜｜λ３｜とする。
また、面状構造は図５に示すように、３つのうち１つの固有値の値が大きく、２つが０に近い特徴を持つことが知られている。例えば、式（１）の固有値は、面状構造からなる対象組織に対して、式（７）のような関係を有する。
また、点状構造は、３つの固有値の全てが大きい特徴を持つことが知られている。例えば、式（１）の固有値は、点状構造からなる対象組織に対して、式（８）のような関係を有する。
したがって、固有値から線状構造らしさ、面状構造および点状構造らしさを判別することができ、判別結果を用いて、３次元画像Ｍ０において、線状構造物である気管支領域、および面状構造物である骨領域を分割することができる。
中空球モデル関数ｆ（ｒ）について、図３Ａ、３Ｂを用いて、以下説明する。図３Ａは、Ｘ方向における中空球モデル関数ｆ（ｒ）のレスポンスを実線で示し、中空球モデル関数ｆ（ｒ）のｘ方向の２階偏微分のレスポンスを破線で示している。図３Ｂの中央部分は、中空球モデル関数ｆ（ｒ）のレスポンスを、ｚ軸方向に等間隔に切断したｘｙ平面図により示している。また、図３Ｂの右部分は、中空球モデル関数ｆ（ｒ）のｘ方向の２階偏微分のレスポンスをｚ軸方向に等間隔に切断したｘｙ平面図により示している。なお、図３Ｂ中央部分、右部分とも各レスポンスを表す平面図をｚ座標値の大きい順に垂直方向上から並べて示している。図３Ｂのｘｙ平面図において、明度が大きい（白い）ほどフィルタ係数が正の大きい値となり、明度が小さい（黒い）ほどフィルタ係数が負の小さい値（フィルタ係数の絶対値の大きい値）となる。
中空の構造物に対して図３Ａに示すような中空球モデル関数を適用してフィルタリングを行った場合には、中空の構造物の輪郭部分が中空球モデルの球壁（半径ｒ＝−Ｒ，Ｒとなる部分）とそれぞれ一致した場合、式（２）により得られるレスポンスが好適に大きくなるため、非常に精度よく判別を行うことができる。このため、検出対象構造物を横断する線分、例えば気管支などの管状構造物の直径などに対して半径Ｒが一致するように、３次元多重解像度画像Ｍｓｉの解像度を設定することが好ましい。
評価部３４は、フィルタリング部３２が算出したヘッセ行列を固有値分解し、３つの固有値λ１，λ２，λ３（ただし、｜λ１｜≦｜λ２｜≦｜λ３｜とする）を算出する。そして、下記の式（１０）、（１１）に示すように、３次元多重解像度画像Ｍｓｉの各画素における線状構造らしさの評価値Ｌ０（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領域分割方法を用いて、３次元画像Ｍ０を対象領域（気管支領域）と背景領域とに分割する（ステップＳＴ５）。そして、表示部５０が分割された対象領域および背景領域をボリュームレンダリング表示し（ステップＳＴ６）、処理を終了する。
このように、本実施形態によれば、中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部３２と、算出された評価行列に基づいて、各画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出する評価部３４とを備えたため、判別対象物の構造的特徴を好適に表す中空球を表す関数により画像のフィルタリングを行うことができる。このため、好適なフィルタリングの結果得られた評価行列から算出される評価値に基づいて、画像に含まれる中空の構造物の点状構造、線状構造および面状構造のうち評価されたものの評価精度を向上することができる。
また、上記第１の実施形態においては、複数サイズの半径Ｒの中空球を表す関数について、該各中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列をそれぞれ算出し、この複数の評価行列に基づいて、それぞれ各画素位置における線状構造らしさおよび面状構造らしさの評価値を算出している。このため、各カーネルサイズに対応するサイズの構造物を好適に評価することができる。また、この評価値に基づいて、上記の第１の実施形態のように対象領域を分割した場合には、好適に複数サイズの中空の構造物を対象領域として分割することができる。
なお、異なるサイズの構造物に対してヘッセ解析を行う手法として、本実施形態のように複数の解像度の画像に対して、一定のカーネルサイズσおよび一定の中空球モデル関数の半径Ｒを適用してヘッセ解析を行う手法を適用してもよく、一定の解像度の画像に対して、カーネルサイズσおよび中空球モデル関数の半径Ｒを異ならせてヘッセ解析を行う手法を適用しても同効果が得られる。例えば、中空の構造物の長軸や短軸の長さなど、検出対象構造物を横断する線分の長さを試験的に得られたデータなどから予め取得しておき、検出対象構造物を横断する線分の長さと半径Ｒの大きさが一致するように、３次元多重解像度画像Ｍｓｉの解像度、あるいは、カーネルサイズσおよび半径Ｒを設定することが考えられる。
なお、計算負荷の抑制と計算処理の高速化のためには本実施形態のようにフーリエ空間においてフィルタリングを行うことが好適であるが、実空間において対象画像に中空球モデル関数ｆ（ｒ）の２階偏微分およびガウシアンカーネルを畳み込み積分（コンボリューション）することによりフィルタリング行ってもよい。
また、上記第１の実施形態においては、線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０の双方を算出しているが、線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０のいずれか一方のみを算出するようにしてもよい。
また、上記第１の実施形態においては、ヘッセ行列を固有値分解することにより取得した固有値を用いて線状構造らしさの評価値Ｌ０および面状構造らしさの評価値Ｐ０を算出しているが、固有値を得ることなく、ヘッセ行列の要素をそのまま用いて線状構造らしさの評価値および面状構造らしさの評価値を算出するようにしてもよい。
本発明の第２の実施形態として、方向判別のための画像処理装置を以下に説明する。図７は方向判別のための画像処理装置の機能ブロック図を表す。第２の実施形態における画像処理装置は、第１の実施形態における分割部４０と評価部３４と、表示部６０を備えておらず、代わりに方向判別部３５を備えている。この方向判別部３５による方向判別処理以外の処理である、画像取得部１０による画像取得処理と、検出領域設定部２０による検出領域設定処理と、フィルタリング部３２によるフィルタリング処理については、第１の実施形態と同じであり、各機能ブロックの機能も共通している。以下、第１の実施形態と異なる点を中心に説明し、第１の実施形態と同じ点については説明を省略する。
方向判別部３５は、フィルタリング部３２が算出した評価行列を固有値解析して固有値および固有ベクトルｅ１，ｅ２，ｅ３を算出し、算出された固有ベクトルの向きに基づいて線状構造物の主軸方向または面状構造の法線方向を判別する。詳細には、フィルタリング部３２が得たヘッセ行列の固有値λ１，λ２，λ３および固有値ベクトルｅ１，ｅ２，ｅ３を取得し、式（６）に示す関係を満たしていれば、固有値λ１に対応する固有ベクトルｅ１の指す方向を線状構造の軸方向と判別し、式（７）に示す関係を満たしていれば、固有値λ３に対応する固有ベクトルｅ３の指す方向を面状構造の法線方向と判別する。
図８は本実施形態において行われる処理を示すフローチャートである。まず、画像取得部１０が、Ｘ線ＣＴ装置２により撮像された２次元画像から３次元画像Ｍ０を生成する（ステップＳＴ１１）。次に、検出領域設定部２０が、３次元画像Ｍ０を等方化するとともに多重解像度変換して解像度が異なる複数の３次元多重解像度画像Ｍｓｉを生成する（ステップＳＴ１２）。
次いで、判別部３０のフィルタリング部３２が、３次元多重解像度画像Ｍｓｉ毎に、中空球モデル関数およびガウシアンカーネル関数を用いたフィルタリングを行い、各画素位置におけるヘッセ行列を算出する（ステップＳＴ１３）。次いで、方向判別部３５が、ヘッセ解析の結果に基づいて、各画素の局所的な固有値および固有値ベクトルを取得し、先述のように式（４）および式（５）に基づいて線状構造物の軸方向または面状構造物の法線方向を算出し（ステップＳＴ１４）、処理を終了する。
以上のように、第２の実施形態においても、第１の実施形態と同様に判別対象物の構造的特徴を好適に表す中空球を表す関数により画像のフィルタリングを行うことができ、結果として評価行列から算出される固有値および固有ベクトルの値の精度を向上することができる。このため、固有値および固有ベクトルに基づいて中空の構造物の各位置において、線状構造の軸方向または面状構造の法線方向を良好に判別することができる。上記第２の実施形態における画像処理装置により得られた線状構造の軸方向や面状構造の法線方向は、例えば気管支や大腸など管状構造のＣＰＲ(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 (8)

  1. 中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部と、
    前記算出された評価行列に基づいて、前記画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出する評価部とを備えたことを特徴とする画像処理装置。
  2. 中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部と、
    前記算出された評価行列を固有値解析して固有値および固有ベクトルを算出し、算出された固有ベクトルの向きに基づいて線状構造の主軸方向または面状構造の法線方向を判別する方向判別部とを備えたことを特徴とする画像処理装置。
  3. 前記中空球を表す関数は、下記式（３）により表されるものであることを特徴とする請求項１または２記載の画像処理装置。
    ただし、ｘ、ｙ、ｚは実空間における３次元座標であり、Ｒは中空球の半径である。
  4. 前記フィルタリング部は、複数サイズの前記中空球を表す関数について、該各中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列をそれぞれ算出するものであることを特徴とする請求項１から３のいずれか１項記載の画像処理装置。
  5. 前記フィルタリング部は、前記中空球を表す関数の二階偏微分行列を用いたフィルタリングを、フーリエ空間において実施するものであることを特徴とする請求項１から４のいずれか１項記載の画像処理装置。
  6. 前記画像が医用画像であり、前記中空の構造物が気管支であることを特徴とする請求項１から５のいずれか１項記載の画像処理装置。
  7. 中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出し、
    前記算出された評価行列に基づいて、前記画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出することを特徴とする画像処理方法。
  8. コンピュータを、
    中空の構造物を表す画像中の各画素位置において、中空球を表す関数の二階偏微分行列を用いたフィルタリングを施した評価行列を算出するフィルタリング部と、
    前記算出された評価行列に基づいて、前記画素位置における点状構造らしさ、線状構造らしさおよび面状構造らしさのうち少なくとも１つ以上の評価値を算出する評価部として機能させることを特徴とする画像処理プログラム。
JP2012055599A 2012-03-13 2012-03-13 画像処理装置および方法並びにプログラム Active JP5748693B2 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2012055599A JP5748693B2 (ja) 2012-03-13 2012-03-13 画像処理装置および方法並びにプログラム
PCT/JP2013/001528 WO2013136750A1 (ja) 2012-03-13 2013-03-08 画像処理装置および方法並びにプログラム
US14/481,741 US9307948B2 (en) 2012-03-13 2014-09-09 Image processing apparatus, method, and program

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2012055599A JP5748693B2 (ja) 2012-03-13 2012-03-13 画像処理装置および方法並びにプログラム

Publications (2)

Publication Number Publication Date
JP2013188289A JP2013188289A (ja) 2013-09-26
JP5748693B2 true JP5748693B2 (ja) 2015-07-15

Family

ID=49160682

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2012055599A Active JP5748693B2 (ja) 2012-03-13 2012-03-13 画像処理装置および方法並びにプログラム

Country Status (3)

Country Link
US (1) US9307948B2 (ja)
JP (1) JP5748693B2 (ja)
WO (1) WO2013136750A1 (ja)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6334902B2 (ja) 2012-11-30 2018-05-30 キヤノンメディカルシステムズ株式会社 医用画像処理装置
EP3563340A2 (en) * 2017-03-09 2019-11-06 St. Jude Medical International Holding S.à r.l. Detection of fiducials in a clinical image

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7020314B1 (en) * 2001-11-13 2006-03-28 Koninklijke Philips Electronics N.V. Black blood angiography method and apparatus
AU2005329247A1 (en) * 2005-03-17 2006-09-21 Algotec Systems Ltd. Bone segmentation
US20110064289A1 (en) * 2009-09-14 2011-03-17 Siemens Medical Solutions Usa, Inc. Systems and Methods for Multilevel Nodule Attachment Classification in 3D CT Lung Images
US8724866B2 (en) 2009-09-14 2014-05-13 Siemens Medical Solutions Usa, Inc. Multi-level contextual learning of data
JP2011104206A (ja) * 2009-11-19 2011-06-02 Hiroshima City Univ 動脈瘤候補の検出支援装置および検出方法
JP5037705B2 (ja) * 2010-03-11 2012-10-03 富士フイルム株式会社 画像処理装置および方法並びにプログラム
JP5265810B2 (ja) * 2010-10-08 2013-08-14 パナソニック株式会社 超音波診断装置、及び体内観察方法
JP5701138B2 (ja) * 2011-04-19 2015-04-15 富士フイルム株式会社 医用画像処理装置および方法、並びにプログラム

Also Published As

Publication number Publication date
US20140376796A1 (en) 2014-12-25
WO2013136750A1 (ja) 2013-09-19
JP2013188289A (ja) 2013-09-26
US9307948B2 (en) 2016-04-12

Similar Documents

Publication Publication Date Title
Sun et al. Automated 3-D segmentation of lungs with lung cancer in CT data using a novel robust active shape model approach
US9183637B2 (en) Image processing apparatus, method, and program
US7684602B2 (en) Method and system for local visualization for tubular structures
JP6570145B2 (ja) 画像を処理する方法、プログラム、代替的な投影を構築する方法および装置
EP2372646B1 (en) Image processing device, method and program
US9563978B2 (en) Image generation apparatus, method, and medium with image generation program recorded thereon
US20170011534A1 (en) Generating a synthetic two-dimensional mammogram
US20050196024A1 (en) Method of lung lobe segmentation and computer system
EP2856428B1 (en) Segmentation highlighter
Bhateja et al. An improved medical image fusion approach using PCA and complex wavelets
JP6458166B2 (ja) 医用画像処理方法及び装置及びシステム及びプログラム
Besler et al. Bone and joint enhancement filtering: Application to proximal femur segmentation from uncalibrated computed tomography datasets
JP5748693B2 (ja) 画像処理装置および方法並びにプログラム
EP1447772B1 (en) A method of lung lobe segmentation and computer system
WO2015111308A1 (ja) ３次元医用画像表示制御装置およびその作動方法並びに３次元医用画像表示制御プログラム
WO2006055031A2 (en) Method and system for local visualization for tubular structures
Shah et al. Quantification and visualization of mri cartilage of the knee: A simplified approach
JP2010200925A (ja) 画像処理装置および方法並びにプログラム
Gill et al. Segmentation of lungs with interstitial lung disease in CT Scans: A TV-L 1 based texture analysis approach
Zhou et al. Vessel boundary extraction using ridge scan-conversion deformable model
JP2014161388A (ja) 画像処理装置、画像処理方法、画像処理装置の制御プログラム、記録媒体
Blons et al. PerceptFlow: Real-Time Ultrafast Doppler Image Enhancement Using Deep Convolutional Neural Network and Perceptual Loss
Ma et al. 3D Reconstruction of Carotid Artery from Ultrasound Images
CN116758158A (zh) 一种医学超声检测探头位姿调整方法、介质及系统
Rodrigues et al. A region-based algorithm for automatic bone segmentation in volumetric CT

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20140526

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20150421

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20150512

R150 Certificate of patent or registration of utility model

Ref document number: 5748693

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250