JP5879291B2

JP5879291B2 - 画像処理装置、画像処理プログラムおよび画像処理装置の作動方法

Info

Publication number: JP5879291B2
Application number: JP2013066092A
Authority: JP
Inventors: 嘉郎北村
Original assignee: Fujifilm Corp
Current assignee: Fujifilm Corp
Priority date: 2013-03-27
Filing date: 2013-03-27
Publication date: 2016-03-08
Anticipated expiration: 2033-03-27
Also published as: DE112014001697T5; US9965698B2; WO2014156089A1; US20160019435A1; JP2014191564A

Description

本発明は、グラフカットやＱＰＢＯ（Quadratic Pseudo-Boolean Optimization）などの離散最適化手法を用いて複数の領域に画像を分割する画像処理装置、画像処理プログラムおよび画像処理装置の作動方法に関する。

近年、グラフの最小切断（グラフカット）アルゴリズムを使ったエネルギー最小化が、画像処理に盛んに応用されるようになってきた。特に画像のセグメンテーション（領域分割）の問題をエネルギー最小化問題として効率的に解く方法が提案されている。

グラフカットは、２変数に依存する２次のエネルギーを基本として発展してきたが、Kolmogorovらは３変数に依存する劣モジュラな３次のエネルギーを最小化する方法を提案している（例えば、非特許文献１）。

さらに、４次以上の高次のエネルギー関数を２次のエネルギー関数に変換して解く方法が提案されている。変換後のエネルギー関数の２次項の係数が全て負の場合をサブモジュラと呼び、Max-Flow/Min-Cutアルゴリズムにより高速に最小解を得ることができる。一方、変換後のエネルギー関数の２次項の係数に正の値が含まれる場合（非サブモジュラ）は、ＱＰＢＯアルゴリズムを用いて最小化または近似解を得ることで、元の高次エネルギー関数の最適化が可能である。このエネルギー関数の変換は複数知られており、例えばRotherらによりこの変換が提案されている（例えば、非特許文献２）。

また、グラフカットの手法を用いて、抽出対象となる対象物の領域を抽出するために、対象物の輪郭を表す閉曲線との距離に基づいてエネルギーを設定し、オブジェクト領域とそれ以外の領域とにセグメンテーションする手法を提案したものがある（例えば、特許文献１）。

V. Kolmogorov and R. Zabih, What energy functions can be minimized via Graph Cuts?, IEEE PAMI, 26(2), 2004. C. Rother, P. Kohli, W. Feng, j. Jia, Minimizing Sparse Higher Order Energy Functions of Discrete Variables, CVPR, 2009.

特開２０１２-０２７７１３号公報

非特許文献１および非特許文献２には、高次エネルギーを２次のエネルギーに変換して解を得る方法が示され、特許文献１には、輪郭に基づいてエネルギーを設定することで特定の輪郭を持つ対象物領域とそれ以外の領域とにセグメンテーションする手法が示されている。非特許文献２には、高次のエネルギーを利用してセグメンテーションする一つの応用例として、矩形のウィンドウ内の画素に対して高次のエネルギーを設定している。この方法では高次エネルギーを十分に疎にすることが難しく、計算量が大きいという問題があった。

そこで、本発明では、セグメンテーションに有効な高次のエネルギーを利用した画像処理装置、画像処理プログラムおよび画像処理装置の作動方法を提供することを目的とする。

本発明の画像処理装置は、エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第１の変数または前記対象物領域を除く非対象領域を表す第２の変数にラベル付けをする画像処理装置であって、
対象物領域の輪郭に相似した形状を有する、対象物領域側に存在する第１輪郭、および、前記対象物領域の輪郭を挟んで第１輪郭と対向した位置に存在する非対象物領域側に存在する第２輪郭、を特定する輪郭特定手段と、第１輪郭の全体を形成するｉ(≧1)個の画素と、第２輪郭の全体を形成する（Ｎ−ｉ）(Ｎ≧４、Ｎ−ｉ≧1)個の画素を選択する画素選択手段と、第１輪郭の前記ｉ個の画素が全て対象物領域に属し、かつ、第２輪郭の（Ｎ−ｉ）個の画素が全て非対象領域に属するという条件を満たす場合のＮ次のエネルギーを、前記条件を満たさない場合のＮ次のエネルギーよりも小さくなるように設定するエネルギー設定手段と、Ｎ次のエネルギーを最小化してラベル付けをするラベリング手段とを備えることを特徴とするものである。

本発明の画像処理プログラムは、コンピュータを、エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第１の変数または前記対象物領域を除く非対象領域を表す第２の変数にラベル付けをする画像処理装置として機能させるためのプログラムであって、
対象物領域の輪郭に相似した形状を有する、対象物領域側に存在する第１輪郭、および、前記対象物領域の輪郭を挟んで第１輪郭と対向した位置に存在する非対象物領域側に存在する第２輪郭、を特定する輪郭特定手段と、第１輪郭の全体を形成するｉ(≧1)個の画素と、第２輪郭の全体を形成する（Ｎ−ｉ）(Ｎ≧４、Ｎ−ｉ≧1)個の画素を選択する画素選択手段と、第１輪郭のｉ個の画素が全て対象物領域に属し、かつ、第２輪郭の（Ｎ−ｉ）個の画素が全て非対象領域に属するという条件を満たす場合のＮ次のエネルギーを、条件を満たさない場合のＮ次のエネルギーよりも小さくなるように設定するエネルギー設定手段と、Ｎ次のエネルギーを最小化してラベル付けをするラベリング手段として機能させるものであることを特徴とするものである。

本発明の作動方法は、輪郭特定手段と、画素選択手段と、エネルギー設定手段と、ラベリング手段とを有し、エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第１の変数または前記対象物領域を除く非対象領域を表す第２の変数にラベル付けをする画像処理装置の作動方法であって、
輪郭特定手段により、対象物領域の輪郭に相似した形状を有する、対象物領域側に存在する第１輪郭、および、対象物領域の輪郭を挟んで第１輪郭と対向した位置に存在する前記非対象物領域側に存在する第２輪郭、を特定する輪郭特定ステップと、画素選択手段により第１輪郭の全体を形成するｉ(≧1)個の画素と、第２輪郭の全体を形成する（Ｎ−ｉ）(Ｎ≧４、Ｎ−ｉ≧1)個の画素を選択する画素選択ステップと、エネルギー設定手段により第１輪郭の前記ｉ個の画素が全て対象物領域に属し、かつ、第２輪郭の（Ｎ−ｉ）個の画素が全て非対象領域に属するという条件を満たす場合のＮ次のエネルギーを、前記条件を満たさない場合のＮ次のエネルギーよりも小さくなるように設定するエネルギー設定ステップと、ラベリング手段により、Ｎ次のエネルギーを最小化してラベル付けをするラベリングステップを備えたことを特徴とするものである。

「対象物領域の輪郭に相似した形状」とは、対象物領域の輪郭の形状に似通った形状を有するものをいい、第１輪郭は対象物領域の輪郭から略一定の距離分離れて対象物領域内側に位置するものが好ましく、第２輪郭は対象物領域の輪郭から略一定の距離分離れて対象物領域外側（非対象領域側）に位置するものが好ましい。

また、輪郭特定手段が、第１輪郭と第２輪郭とが所定の距離離間するように第１輪郭および第２輪郭を特定するものが好ましい。

また、画像データを用いて前記対象物領域の輪郭を推定する輪郭推定手段をさらに備え、輪郭特定手段が、推定した輪郭に基づいて、前記第１輪郭および前記第２輪郭を特定するものであってもよい。

また、輪郭推定手段が、エッジを検出可能な微分フィルタリングを用いて画像データから対象物領域の輪郭を推定するものであってもよい。

また、輪郭推定手段が、所定のパラメータで表現された形状モデルを用いて前記画像データから前記対象領域の輪郭を推定するものであってもよい。

さらに、エネルギー設定手段が、条件を満たさない場合のＮ次のエネルギーを全て等しい値に設定するものが望ましい。

さらにまた、ラベリング手段が、エネルギーの最小化のためにＱＰＢＯアルゴリズムを用いるものであってもよい。

第１の画像処理装置の構成を示したブロック図画像処理の流れを示すフローチャート心臓の造影CT画像から三次元のCPR画像を生成する手法を説明するための図心臓の造影CT画像から三次元のCPR画像を生成した例２クラスに分類する場合のグラフの一例高次のエネルギーを擬ブール関数で変換した２次式を表すグラフの一例血管の断面形状を説明するための図血管領域が画像上にあらわれる形状の例血管領域が背景領域から分離されるパターンを説明するための図２次のエネルギーの設定を説明するための図対象物領域の輪郭と第１輪郭と第２輪郭の関係を示す図第１輪郭と第２輪郭を離間させた一例第１輪郭上の画素と第２輪郭上の画素を選択する手法を説明するための図三次元OOFフィルタを用いて検出した血管らしい円とＮ次のエネルギーの設定方法を説明するための図三次元CPR画像中から検出された複数の円から血管領域をセグメンテーションする手法を説明するための図従来の手法による血管領域の抽出と本発明の手法による血管領域の抽出の違いを説明するための図Ｎ次のエネルギーを設定する形状とその形状を構成する画素のとり方を説明するための図第２の画像処理装置の構成を示したブロック図対象領域の輪郭の円を複数の円弧で定義する手法を説明するための図

以下、図面を参照して本発明の実施形態を詳細に説明する。本発明の実施形態では、グラフカットの手法を用いて画像中の各画素を対象物領域を表す第１の変数または対象物領域を除く非対象領域を表す第２の変数に領域分割する手法について説明する。本発明の画像処理装置は、画像処理プログラムがコンピュータにロードされて実行されることにより実現する。また、画像処理プログラムはＣＤ−ＲＯＭ等の記憶媒体に記憶されて配布され、ＣＤ−ＲＯＭ等の記憶媒体からコンピュータにインストールされる。あるいは、インターネット等のネットワークを介してプログラムが配布されて、コンピュータにインストールされる。

図１に、本発明に係る第１の画像処理装置１の構成を示したブロック図を示し、図２に、本発明に係る画像処理装置１の処理の流れを示したフローチャートを示す。

図１に示すように、画像処理装置１は、画像データ入力受付手段１１と、前処理手段１２と、輪郭特定手段１３と、画素選択手段１４と、エネルギー設定手段１５と、ラベリング手段１６と、表示手段１７とを備えている。

画像データ入力受付手段１１は、画像処理対象の画像データＰの入力を受け付け、記憶装置に記憶する（S1）。画像データＰは、単純Ｘ線撮影画像や、ＣＴ装置やＭＲＩ装置で撮影された断層画像等である。また、画像データＰには、複数の臓器や組織あるいは病変部が撮影されている。

以下、画像データＰが心臓の造影CT画像（ボクセルデータ）であり、この心臓の造影CT画像から冠動脈の輪郭を検出する処理について説明する。

前処理手段１２は、既知の冠動脈抽出処理により、図３に示すように、造影CT画像から冠動脈の中心線ｌ（破線）を抽出する（S2）。冠動脈抽出処理は、特許第４７０９２９０号や特許第４７１７９３５号に示される方法をはじめ、種々の方法が提案されている。造影CT画像から血管の中心線に直行する断面画像Ｓ_１，Ｓ_２，Ｓ_３を生成し（図３の左図）、この断面画像の血管の中心座標が直線Ｌ上に位置するように断面画像Ｓ_１，Ｓ_２，Ｓ_３を配置して（図３の右図）、三次元のCPR画像を生成する（S3）。図４に心臓の造影CT画像（同図（Ａ）、白い血管に対して処理を行う）から三次元のCPR画像（同図（Ｂ）（Ｃ）、ただし（Ｃ）は血管の断面の図）を生成した例を示す。

前処理手段１２で生成した三次元のCPR画像を背景領域、血管領域の２つに領域にセグメンテーションする。

まず、画像を２つに領域に分割するために用いるグラフカットの手法について説明する。画像のセグメンテーションは、画像の各画素に対応する変数｛ｘ_１，ｘ_２，・・・，ｘ_ｎ｝、ｘ∈｛１，０｝と定義したとき、次式（１）で表される１次のエネルギーＥ^ｉと２次のエネルギーＥ^ｉｊからなるエネルギー関数を最小化することで解くことができる。
Ｖは画像を構成する画素の集合であり、Ｎiは画素iと隣接関係がある画素の集合である。

式（１）の１次のエネルギーＥⁱ(ｘi)は、各画素に割り当てられるラベルのみに依存する値であり、各画素にどのラベルを割り当てたかによって決まる直接的な影響が現れる。第２項目の和は、隣接する画素に与えられるラベルがどのような関係にあるべきかという事前知識を反映するようにエネルギーを定義したものである。

図５に示すように、式（１）を表した有向グラフを設定する。グラフは各画素に対応する頂点ｘ_i（または、頂点ｘ_ｊ）と、さらに各ラベルに対応する２つの頂点ｓ（＝０）と頂点ｔ（＝１）を持ち、頂点ｘ_i、頂点ｘ_ｊそれぞれについて、頂点ｓから頂点ｘ_iへ向かうエッジと頂点ｘ_iから頂点ｔへ向かうエッジと、隣接する画素に対応する頂点ｘ_ｉ、ｘ_ｊ間のエッジで構成される。さらに、頂点ｓから頂点ｘ_iへ向かうエッジと頂点ｘ_iから頂点ｔへ向かうエッジに対して、１次のエネルギーＥ^ｉ(ｘ_ｉ)を定義し、隣接する画素に対応する頂点を結ぶエッジに対しては、２次のエネルギーＥ^ｉｊ(ｘ_ｉ,ｘ_ｊ)を定義する。式（１）であらわされるエネルギーＥ(ｘ)の最小化により各画素に対するラベルを得るためには、最小エネルギーを与えるラベルと最小切断が対応するように各エッジのエネルギー（重み）を定めればよい。

例えば、この２次のエネルギーＥ^ｉｊ(ｘ_ｉ,ｘ_ｊ)が下式で与えられるとする。

このとき２次のエネルギーＥ^ｉｊを擬ブール式で書くと、下のような２次式となる。
擬ブール式において、２次の項の係数が全て負の場合をサブモジュラと呼び、Max-Flow/Min-Cutアルゴリズムにより高速に最小解を得ることができる。一方、２次の項の係数に正の値が含まれる場合（つまり、非サブモジュラ）には、ＱＰＢＯアルゴリズムを用いて最小化または準最適解を得ることができる（詳細は、C. Rother, V. Kolmogorov, V. Lempitsky, M. Szummer, Optimizing Binary MRFs via Extended Roof Duality, CVPR, 2007.を参照）。

さらに、３次以上のエネルギー関数については、３次以上のエネルギー関数を２次のエネルギー関数に変換した上で、変換後の２次のエネルギー関数を従来から知られているＱＰＢＯアルゴリズムで最適化することにより、元の３次以上の高次エネルギー関数の最適化が可能である。

例えば、Rotherらは補助変数ｚ_０，ｚ_１を追加してグラフを変換する手法を開示している。例えば、３次のエネルギー関数が、下式で与えられているとする。

この３次のエネルギーＥ^ｉｊｋを擬ブール式で書くと、下のような２次式となる。
この変換をグラフで表すと、図６のようになる。このように変換した後、補助変数を含めて最小化することで解を得る。その結果得られる解は、もとの高次関数の解と等しい（詳細は、C. Rother, P. Kohli, W. Feng, j. Jia, Minimizing Sparse Higher Order Energy Functions of Discrete Variables, CVPR, 2009.を参照）。

４次以上のエネルギーについても同様の変換を行うことで解を得ることが可能である。この変換は、高次エネルギーの組合せ表のうち値が異なる要素ごとに、変換を行う必要がある。複雑な場合には追加する補助変数の数が多くなり最適化が難しくなるが、エネルギーが疎の場合、例えば、式（４）のように特定の値のみが小さい場合は、相対的にエネルギーが小さい要素だけを変換すればよい。従って変換後の問題が比較的単純で、実用的な計算時間で解が得られる。

以下、４次以上のエネルギーを用いて、形状が既知である対象物が存在する対象物領域とそれ以外の領域（非対象物領域）にセグメンテーションする手法について説明する。以下、対象物を血管とし、非対象物を背景とする。

図７に示すように、血管の断面形状は概ね円形であり、画像上の血管の輪郭に円（破線）をあてはめたとき、その円の輪郭の内側が血管領域であり、外側が背景領域である。

２次のエネルギーで表現する場合は、２つの画素のうち背景領域に位置する画素が０をとり、血管領域に位置する画素が１を取る組合せのみ、小さなエネルギーを設定してグラフカットを行うことで対象物領域と背景領域に分離することができる。

図８に示すように、エネルギーを２つの円の輪郭が互いに重なった状態で設定する場合を考える。このとき、円の輪郭に２次のエネルギーを設定して血管領域と背景領域を分離すると、図９に示す６通りのパターンに分離される可能性がある。図９（Ａ）のように、左または右の円のどちらかの内側（白い部分）を血管として背景（黒い部分）から分離される場合と、図９（Ｂ）のように円が一部欠けた形状のどちらかの内側（白い部分）を血管として分離される場合と、図９（Ｃ）に示す形状で分離される場合がある。２次のエネルギーは、２つの画素が円の輪郭を跨ぐ（図１０参照）場合に小さいエネルギーを設定したとする。図９（Ａ）（Ｂ）のいずれのパターンでも、エネルギーの和は円の輪郭の長さに比例して小さくなるので、輪郭の長さが等しい図９（Ａ）（Ｂ）のいずれのパターンの分離も同じ確率で起こる可能性がある。さらに、図９（Ｃ）のようなパターンでもエネルギーの和は相対的に小さいので、解の候補の一つである。しかし、血管は円に近い形状をもつことを考慮すると、図９（Ｂ）（Ｃ）のようなパターンで分離されるのは好ましくない。

そこで、円に近い形状で血管領域をセグメンテーションするために、１次のエネルギーと２次のエネルギーに加えて、Ｎ（≧４）個の画素の画素値を変数とするＮ次のエネルギーを用いる。Ｎ次のエネルギーには、血管に近い大きさを持った円形の輪郭の周囲の背景に位置する画素が０（図７の●）をとり、輪郭の周囲の血管に位置する画素が１（図７の○）をとる組合せのみ小さいエネルギーを設定する。円は互いに重なりあう部分があるので、２つの円のＮ次エネルギーは一方が小さくなるともう一方が小さくなることはない。このように排他的に一方の円が選ばれるように作用するので、血管領域が円に近い形状になるように画像がセグメンテーションされる。

４次以上の高次エネルギーを加えたエネルギー関数Ｅを次式のように定義する。
ここでｃは高次エネルギーを構成する画素の集合であり、Ｘｃは関係する画素のバイナリ変数（０または１）からなるベクトルである。Ｃはｃの集合を表す。なお、各ｃの次数はそれぞれ異なってもよく、例えば大きな円に関連する画素の数は多く、小さな円に対しては少ない。

まず、輪郭特定手段１３で、血管（対象物）領域の内側の輪郭と外側の輪郭を特定する。具体的には、図１１に示すように、血管の輪郭（実線）を挟んだ血管（対象物）領域側に存在する第１輪郭（破線）と、血管の輪郭を挟んで第１輪郭と対向した位置に存在する背景（非対象物）領域側に存在する第２輪郭（一点鎖線）を特定する。第１輪郭と第２輪郭は、血管の輪郭に相似した形状にする。血管の輪郭は略円形であるので、血管の大きさに近い円を設定し、さらに、その内側に第１輪郭の円を設け、外側に第２輪郭の円を設ける。

しかし、完全に血管の大きさをもった円形に固定してセグメンテーションが行われると、柔軟性がなく都合が悪い。実際の画像にあらわれる血管は、完全な円形ではなく少し歪んだ円になることもあり、実際の画像に即した輪郭線が得られるようにしたい。そこで、図１２に示すように、血管領域側の第１輪郭（内側の円）と背景領域側の第２輪郭（外側の円）との間を所定の距離離して隙間を設ける（S4）。こうすることで、その間であれば自由な輪郭（破線）をとっても設定した高次エネルギーの条件は成立し、エネルギーが小さくなる。すなわち高次が成立する範囲で、２次のエネルギーが最小となるように柔軟に輪郭線が選択される。また、完全な円を設定するのではなく、図１９に示すように、第１輪郭と第２輪郭の円を分割して複数の円弧を設定するようにしてもよい。このように分割すると、円の一部のパターンが一致するだけでもエネルギーを小さくする効果があるため、より柔軟なセグメンテーション結果を得ることができる。

次に、画素選択手段１４は、三次元CPR画像から血管らしい形状が現れている場所を探して、図１３に示すように、第１輪郭上のｉ個の画素と第２輪郭上のＮ-ｉ個の画素を選択する。

エネルギー設定手段１５は、グラフに１次のエネルギーと２次のエネルギーとともに高次のエネルギーを設定する。１次のエネルギーは頂点ｓ（＝０：背景）から頂点ｘ_iへ向かうエッジと頂点ｘ_iから頂点ｔ（＝１：血管）へ向かうエッジとに設定され、各画素の値が血管らしい値を持っている場合は、頂点ｘ_iから頂点ｔへ向かうエッジより頂点ｓから頂点ｘ_iへ向かエッジに小さいエネルギーを設定し、各画素の値が背景らしい値を持っている場合は、頂点ｓから頂点ｘ_iへ向かうエッジより頂点ｘ_iから頂点ｔへ向かうエッジに小さいエネルギーを設定する（S5）。２次のエネルギーは隣接する画素に対応する頂点ｘ_ｉ、ｘ_ｊ間のエッジに設定され、隣接画素値の差分が大きいほど小さいエネルギーを設定する（S6）。

さらに、エネルギー設定手段１５は、高次のエネルギーを効率的に設定するために、三次元CPR画像中の円らしい形状があらわれている位置を特定して高次のエネルギーを設定する。第１輪郭のｉ個の画素が全て血管らしい画素値を持ち、かつ、第２輪郭の（Ｎ−ｉ）個の画素が全て背景らしい画素値を持つ場合は輪郭としての条件を満たすので、この条件を満たす場合は、高次のエネルギーは小さいエネルギー（その重みは係数が負で絶対値が大きい）を設定し、この条件を満たさない場合には、Ｎ次のエネルギーはこの条件を満たした場合のエネルギーより大きくなるように設定する。

円形のＮ次のエネルギーを設定する位置を特定するために、三次元OOFフィルタを用いる。三次元OOFフィルタは半径Ｒの球状（２次元断面では円形）のフィルタ形状を持ち、図１４（Ａ）に示すように、半径Ｒの円（破線）の形状が現れているところにレスポンスがあらわれる。そこで、三次元CPR画像Ｐの各位置に三次元OOFフィルタを施して、三次元CPR画像から円形らしい形状があらわれている位置を特定する。また、三次元OOFフィルタは球の表面における２階偏微分を算出し、その２階偏微分行列を固有値分解して得た３つの固有値から、円柱構造物らしさを算出することができる。三次元OOFフィルタを施した位置に円柱構造物が存在すれば、３つの固有値のうち２つが大きく、１つがゼロに近い値をとる。同時に得られる３つ固有ベクトルのうち、最もゼロに近い固有値を持つ固有ベクトルが円柱構造物の主軸ｔ（走行方向）に相当する。この主軸と直交する平面Ｑ上に、所定の大きさを持つ血管の輪郭らしい画像が存在している（三次元OOFフィルタの詳細は、M. Law et. al., “Three Dimensional Curvilinear Structure Detection Using Optimally Oriented Flux”, Proceedings of ECCV, pp368-382, 2008.を参照）。フィルタリング処理によって、このような円が現れている平面Ｑを検出する（S7）。また血管の正確な大きさは事前にわからないので、様々な大きさの円に対応したOOFフィルタで検出を行う。図１５に三次元CPR画像の各位置に三次元OOFフィルタを施して複数の位置で円が検出された様子を示す。

図１４（Ｂ）に示すように、画素選択手段１４はこの平面Ｑ上でＮ個の画素を選択して、エネルギー設定手段１５でＮ次のエネルギーを設定する（S8）。三次元OOFフィルタで三次元CPR画像Ｐをスキャンし、そのレスポンスが高い位置の平面Ｑにおいて、選択したＮ個の画素が輪郭の条件を満たしたとき相対的に小さい値になるように高次のエネルギーを設定する。一方、三次元OOFフィルタのレスポンスが低い位置の平面Ｑで選択したＮ個の画素が輪郭の条件を満たす場合には、レスポンスが高い位置の平面Ｑで選択したＮ個の画素が輪郭の条件を満たす場合より大きくなるようにＮ次のエネルギーを設定し、かつ、円の条件を満たしていないときの高次のエネルギーよりは小さいエネルギーを設定する。例えば、三次元OOFフィルタのレスポンスが高い位置では、輪郭の条件を満足したときのＮ次のエネルギーに大きい負のエネルギー（係数が負で絶対値が大きい）を設定し、三次元OOFフィルタのレスポンスが低い位置では、輪郭の条件を満足したときのＮ次のエネルギーに小さい負のエネルギー（係数が負で絶対値が小さい）を設定する。

以上のように適応的にエネルギーの大きさを変える理由は、フィルタのレスポンスが高いほど円が存在する可能性が高いと事前に知っているからである。一方で、データから統計的に算出することも可能である。例えば、血管がどのようにセグメンテーションされるべきかを示す正解のデータがある場合には、三次元OOFフィルタのレスポンスに対してＮ次のエネルギーが小さくなるケースの確率を統計的に算出（学習）しておき、その確率に応じてエネルギーの重みを設定するようにしてもよい。この確率は、上述の処理によって設定するＮ次のエネルギーが小さくなるラベルの組み合わせが、正解のデータのラベルと一致する割合として算出できる。

また、エネルギー設定手段１５が、輪郭の条件を満たさない場合のＮ次のエネルギーを全て等しい値に設定するようにしてもよい。このように輪郭の条件を満たさない場合のエネルギーを等しい値にすれば、エネルギーが疎になるので少ない付加変数で高次を低次に変換することができ、解が比較的高速に得られる。

ラベリング手段１６は、設定された高次（３次以上）のエネルギーを２次のエネルギー関数に変換した上で、変換後の２次のエネルギー関数を従来から知られているＱＰＢＯアルゴリズムで最適化して、血管領域と背景領域の２つの領域にセグメンテーションする（S9）。

ここで、円形の高次エネルギーを設定した時の効果について説明する。図１６（Ａ）の２次元の原画像に対して、従来の２次元のエネルギーを設定してグラフカット処理を行った場合には、図１６（Ｂ）に示すように、かなり歪んだ円形で抽出される。一方、図１６（Ｃ）に示すような円形の高次エネルギーを設定した場合には、図１６（Ｄ）に示すような円形に近い形状にセグメンテーションが行われる。

このように円形の高次エネルギーを設定することにより、図１５に示すように抽出された円形（図１５（Ａ））が、３次元的に連続するよう血管領域がセグメンテーションされる（図１５（Ｂ））。また円形から逸脱した形状でセグメンテーションされることを防ぐことができる。

表示手段１７は、ラベリング手段１６の結果に応じて、血管領域のみの画像を表示装置上に表示する(S10)。

上述では、血管のように輪郭の形状が円形であるものに高次のエネルギーを設定する場合について説明したが、高次エネルギーとして設定する輪郭の形はどんな形でも良い。様々な輪郭を挟んで対象物領域側の第１輪郭と非対象物領域側の第２輪郭上の画素を選択して、第１輪郭の画素が全て対象物領域に属し、かつ、第２輪郭の画素が全て非対象領域に属するという場合の高次のエネルギーを、この条件を満たさない場合の高次のエネルギーよりも小さくなるように設定する。

第１輪郭は、対象物領域の輪郭と相似になるように対象物領域の輪郭から常に略一定の距離内側（対象部領域側）に離れたものが好ましく、第２輪郭は、対象物領域の輪郭と相似になるように対象物領域の輪郭から常に略一定の距離外側（非対象部領域側）に離れたものが好ましいが、予め決められた誤差の範囲内で対象物領域の輪郭から距離が異なるように特定されても良い。

また、図１７に示すように、輪郭を構成するためには４個以上の画素を用い、少なくとも第１輪郭または第２輪郭に少なくとも１つ以上の画素があれば、輪郭の形状を特定することができる。図１７（Ａ）に、直線に沿って対象物領域と非対象物領域に分離される例を示し、図１７（Ｂ）に、角のある形状で対象物領域と非対象物領域に分離される例を示している。他にも三角形や矩形など様々な形状を特定することができる。あるいは、面や曲面の形状を特定して、３次元画像中の面や曲面を検出し、分離するようにしてもよい。

また、上述では、血管のように予め形状が分かっている場合について説明したが、図１８に示す第２の画像処理装置１ａのように、画像データＰを用いて対象物領域の輪郭を推定する輪郭推定手段１８を設ける構成にしても良い。

輪郭推定手段１８は、輪郭を推定する手法として、エッジが検出可能な微分フィルタリングを利用することができる。輪郭特定手段１３では、輪郭の近傍で微分フィルタのレスポンスが正から負に変わるので、微分フィルタの正のレスポンスを持つ位置を第１輪郭（対象物領域）とし負のレスポンスを持つ位置を第２輪郭（非対象物領域）として特定し、画素選択手段１４では第１輪郭と第２輪郭上の画素を選択して、エネルギー設定手段１５で正のレスポンスを持つ位置の第１輪郭と負のレスポンスを持つ位置の第２輪郭に組み合わせに対して最小になるようにエネルギーを設定する。以下、前述の通りラベリング手段１６を用いてエネルギー最小化により対象領域と非対象領域にセグメンテーションを行なった後、対象領域のみの画像にするなど適切な画像を生成して表示手段１７で画像を表示装置上に表示する。

あるいは、ある程度形状が分かっている場合には、輪郭推定手段１８は、所定のパラメータで表現された形状モデルを用いて画像データから対象領域の輪郭を推定することもできる。例えば、ＡＳＭ(Active Shape Model)のように、対象物領域の形状を主成分分析や独立成分分析等を行って平均形状と平均形状を変形するためのベクトルを求めておき、各ベクトルのパラメータを変えることで、平均形状から様々な形状に変形する手法が知られている。この形状モデルにおいて、各種パラメータを変化させながら、対応する輪郭位置に高次エネルギーを設定するようにしてもよい。このとき、輪郭が柔軟に変化してもよいように、非対象領域（背景）上にある第１輪郭と対象領域上にある第２輪郭の画素を一定の間隔をあけて選択するようにするとよりよい。

なお、上述の実施例において円を複数の円弧に分割したように、所定のパラメータで表現された形状モデルを複数の線分に分割し、線分ごとにエネルギーを設定するようにしてもよい。

非特許文献２のRotherらの手法では高次のエネルギーを用いることにより、特定の形状をもつ構造物の領域分割の精度を向上させることを可能にしている。一方、本発明は、上記で説明したように、矩形のウィンドウのように固定せず、画素を柔軟に選択する点に特徴がある。これにより、設定する高次エネルギーを疎にすることができ、少ない計算量で高速に解を得ることが可能である。画素を選択する方法は、セグメンテーション対象の輪郭形状に沿って行うことが好ましい。例えば３次元画像中の血管輪郭は円筒形状をしているから、円形となる画素の組み合わせを選択する。また対象の輪郭形状が統計的なモデル（例えばＡＳＭ）として学習されている場合は、そのモデルの輪郭に沿った画素を選択する。このように対象形状の輪郭に沿った画素を選択することが肝要であり、対象構造はなんでもよい。

本実施の形態では、グラフカットの手法をもちいて、エネルギー最小化について説明したが、離散関数を最小化するものであればよい。また、エネルギー最小化の具体的な手法にはグラフカットやＱＰＢＯ、Tree Reweighted Message Passing（詳細は、V. Kolmogorov, Convergent Tree-reweighted Message Passing for Energy Minimization, IEEE PAMI, 28(10), October 2006.を参照）などがある。また、グラフカットを用いて対象物領域と非対象物領域のそれぞれについてラベル付けする方法については、‘Y. Boykov and V. Kolmogorov, An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision, IEEE PAMI, 26(9), 1124-1137, (2004)’に詳細に説明されている。ＱＰＢＯを用いて対象物領域と非対象物領域のそれぞれについてラベル付けする方法は、非特許文献２に詳細に説明されている。

１、１ａ画像処理装置
１１画像データ入力受付手段
１２前処理手段
１３輪郭特定手段
１４画素選択手段
１５エネルギー設定手段
１６ラベリング手段
１７表示手段
１８輪郭推定手段

Claims

エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第１の変数または前記対象物領域を除く非対象領域を表す第２の変数にラベル付けをする画像処理装置であって、
前記対象物領域の輪郭に相似した形状を有する、前記対象物領域側に存在する第１輪郭、および、前記対象物領域の輪郭を挟んで前記第１輪郭と対向した位置に存在する前記非対象物領域側に存在する第２輪郭、を特定する輪郭特定手段と、
前記第１輪郭の全体を形成するｉ(≧１)個の画素と、前記第２輪郭の全体を形成する（Ｎ−ｉ）(Ｎ≧４、Ｎ−ｉ≧１)個の画素を選択する画素選択手段と、
前記第１輪郭の前記ｉ個の画素が全て前記対象物領域に属し、かつ、前記第２輪郭の（Ｎ−ｉ）個の画素が全て前記非対象領域に属するという条件を満たす場合のＮ次のエネルギーを、前記条件を満たさない場合のＮ次のエネルギーよりも小さくなるように設定するエネルギー設定手段と、
前記Ｎ次のエネルギーを最小化してラベル付けをするラベリング手段とを備える画像処理装置。
前記輪郭特定手段が、前記第１輪郭と前記第２輪郭とが所定の距離離間するように第１輪郭および第２輪郭を特定するものであることを特徴とする請求項１記載の画像処理装置。
前記画像データを用いて前記対象物領域の輪郭を推定する輪郭推定手段をさらに備え、
前記輪郭特定手段が、前記推定した輪郭に基づいて、前記第１輪郭および前記第２輪郭を特定することを特徴とする請求項１または２記載の画像処理装置。
前記輪郭推定手段が、エッジを検出可能な微分フィルタリングを用いて前記画像データから対象物領域の輪郭を推定することを特徴とする請求項３記載の画像処理装置。
前記輪郭推定手段が、所定のパラメータで表現された形状モデルを用いて前記画像データから前記対象領域の輪郭を推定することを特徴とする請求項３記載の画像処理装置。
前記エネルギー設定手段が、前記条件を満たさない場合のＮ次のエネルギーを全て等しい値に設定することを特徴とする請求項１〜５いずれか記載の画像処理装置。
前記ラベリング手段が、エネルギーの最小化のためにＱＰＢＯアルゴリズムを用いることを特徴とする請求項１〜６いずれか記載の画像処理装置。
コンピュータを、エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第１の変数または前記対象物領域を除く非対象領域を表す第２の変数にラベル付けをする画像処理装置として機能させるための画像処理プログラムであって、
前記対象物領域の輪郭に相似した形状を有する、前記対象物領域側に存在する第１輪郭、および、前記対象物領域の輪郭を挟んで前記第１輪郭と対向した位置に存在する前記非対象物領域側に存在する第２輪郭、を特定する輪郭特定手段と、
前記第１輪郭の全体を形成するｉ(≧１)個の画素と、前記第２輪郭の全体を形成する（Ｎ−ｉ）(Ｎ≧４、Ｎ−ｉ≧１個の画素を選択する画素選択手段と、
前記第１輪郭の前記ｉ個の画素が全て前記対象物領域に属し、かつ、前記第２輪郭の（Ｎ−ｉ）個の画素が全て前記非対象領域に属するという条件を満たす場合のＮ次のエネルギーを、前記条件を満たさない場合のＮ次のエネルギーよりも小さくなるように設定するエネルギー設定手段と、
前記Ｎ次のエネルギーを最小化してラベル付けをするラベリング手段として機能させるものであることを特徴とする画像処理プログラム。
輪郭特定手段と、画素選択手段と、エネルギー設定手段と、ラベリング手段とを有し、エネルギー最小化手法を用いて、画像データを構成する各画素を対象物領域を表す第１の変数または前記対象物領域を除く非対象領域を表す第２の変数にラベル付けをする画像処理装置の作動方法であって、
前記輪郭特定手段により、前記対象物領域の輪郭に相似した形状を有する、前記対象物領域側に存在する第１輪郭、および、前記対象物領域の輪郭を挟んで前記第１輪郭と対向した位置に存在する前記非対象物領域側に存在する第２輪郭、を特定する輪郭特定ステップと、
前記画素選択手段により、前記第１輪郭の全体を形成するｉ(≧１)個の画素と、前記第２輪郭の全体を形成する（Ｎ−ｉ）(Ｎ≧４、Ｎ−ｉ≧１)個の画素を選択する画素選択ステップと、
前記エネルギー設定手段により、前記第１輪郭の前記ｉ個の画素が全て前記対象物領域に属し、かつ、前記第２輪郭の（Ｎ−ｉ）個の画素が全て前記非対象領域に属するという条件を満たす場合のＮ次のエネルギーを、前記条件を満たさない場合のＮ次のエネルギーよりも小さくなるように設定するエネルギー設定ステップと、
前記ラベリング手段により、前記Ｎ次のエネルギーを最小化してラベル付けをするラベリングステップを備えた画像処理装置の作動方法。