TOP > 国内特許検索 > 画像処理装置、方法、及びプログラム > 明細書

明細書 :画像処理装置、方法、及びプログラム

発行国 日本国特許庁(JP)
公報種別 再公表特許(A1)
発行日 平成29年6月1日(2017.6.1)
発明の名称または考案の名称 画像処理装置、方法、及びプログラム
国際特許分類 G06T   7/00        (2017.01)
G06T   7/60        (2017.01)
FI G06T 7/00 350B
G06T 7/60 150S
国際予備審査の請求 未請求
全頁数 39
出願番号 特願2016-544238 (P2016-544238)
国際出願番号 PCT/JP2015/073277
国際公開番号 WO2016/027840
国際出願日 平成27年8月19日(2015.8.19)
国際公開日 平成28年2月25日(2016.2.25)
優先権出願番号 2014169911
優先日 平成26年8月22日(2014.8.22)
優先権主張国 日本国(JP)
指定国 AP(BW , GH , GM , KE , LR , LS , MW , MZ , NA , RW , SD , SL , ST , SZ , TZ , UG , ZM , ZW) , EA(AM , AZ , BY , KG , KZ , RU , TJ , TM) , EP(AL , AT , BE , BG , CH , CY , CZ , DE , DK , EE , ES , FI , FR , GB , GR , HR , HU , IE , IS , IT , LT , LU , LV , MC , MK , MT , NL , NO , PL , PT , RO , RS , SE , SI , SK , SM , TR) , OA(BF , BJ , CF , CG , CI , CM , GA , GN , GQ , GW , KM , ML , MR , NE , SN , TD , TG) , AE , AG , AL , AM , AO , AT , AU , AZ , BA , BB , BG , BH , BN , BR , BW , BY , BZ , CA , CH , CL , CN , CO , CR , CU , CZ , DE , DK , DM , DO , DZ , EC , EE , EG , ES , FI , GB , GD , GE , GH , GM , GT , HN , HR , HU , ID , IL , IN , IR , IS , JP , KE , KG , KN , KP , KR , KZ , LA , LC , LK , LR , LS , LU , LY , MA , MD , ME , MG , MK , MN , MW , MX , MY , MZ , NA , NG , NI , NO , NZ , OM , PA , PE , PG , PH , PL , PT , QA , RO , RS , RU , RW , SA , SC , SD , SE , SG , SK , SL , SM , ST , SV , SY , TH , TJ , TM , TN , TR , TT , TZ , UA , UG , US
発明者または考案者 【氏名】清水 昭伸
【氏名】斉藤 篤
出願人 【識別番号】504132881
【氏名又は名称】国立大学法人東京農工大学
個別代理人の代理人 【識別番号】100079049、【弁理士】、【氏名又は名称】中島 淳
【識別番号】100084995、【弁理士】、【氏名又は名称】加藤 和詳
【識別番号】100099025、【弁理士】、【氏名又は名称】福田 浩志
審査請求 未請求
テーマコード 5L096
Fターム 5L096AA09
5L096BA06
5L096DA01
5L096EA39
5L096FA19
5L096FA31
5L096JA09
5L096KA04
要約 受付部102が、入力画像を受け付ける。セグメンテーション部114が、予め計算された固有ベクトルを基底とする固有空間であって、固有空間上の点が、特定の物体の形状の統計的変動を表す統計的形状モデルの形状パラメータを示す固有空間において、固有空間上の点が示す形状パラメータが表す特定の物体の形状の尤もらしさと入力画像中の隣接画素間の画素値の差とに応じた値を表す予め定められた目的関数を最適化するように、入力画像が表す被写体の形状を表す形状パラメータを推定し、推定された形状パラメータが表す特定の物体の形状を事前知識として、被写体の領域を、入力画像から抽出する。これによって、計算量の増大を抑制して、被写体の領域を精度よく抽出することができる。
特許請求の範囲 【請求項1】
特定の物体である被写体を表す入力画像から、前記被写体の領域を抽出する画像処理装置であって、
前記入力画像を受け付ける受付手段と、
(A)前記被写体を表す学習用の複数の画像であって、前記被写体の領域が予め求められた前記複数の画像に基づいて予め計算された固有ベクトルを基底とする固有空間であって、かつ、(B)固有空間上の点が、前記特定の物体の形状の統計的変動を表す統計的形状モデルの形状パラメータを示す固有空間において、
前記入力画像に基づいて、前記固有空間上の点が示す前記形状パラメータが表す前記特定の物体の形状の尤もらしさと前記入力画像中の隣接画素間の画素値の差とに応じた値を表す予め定められた目的関数を最適化するように、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、
前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出するセグメンテーション手段と、
を含む画像処理装置。
【請求項2】
前記セグメンテーション手段は、前記固有空間において、
前記入力画像に基づいて、前記目的関数を最適化するように、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定すると同時に、
前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出する請求項1記載の画像処理装置。
【請求項3】
前記セグメンテーション手段は、
最適な形状パラメータを示す点が表す前記特定の物体の形状を含む形状集合を表す前記固有空間上の凸多胞体を繰り返し分割して最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索する探索アルゴリズムに従って、前記目的関数を最適化するように、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、
前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出する請求項1又は請求項2記載の画像処理装置。
【請求項4】
前記セグメンテーション手段は、前記探索アルゴリズムにおいて、前記凸多胞体に含まれる形状集合に対する前記目的関数の下界を計算することにより、最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索して、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出する請求項3記載の画像処理装置。
【請求項5】
前記セグメンテーション手段は、分割後の2つの凸多胞体の体積が対応するように、前記固有空間上の凸多胞体を分割し、かつ、最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索することを繰り返して、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出する請求項4記載の画像処理装置。
【請求項6】
前記セグメンテーション手段は、前記固有空間上にサンプリング点を設定し、前記固有空間上に設定されたサンプリング点から決定される超平面を用いて前記凸多胞体を分割して最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索することを繰り返して、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、
前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出する請求項4又は請求項5記載の画像処理装置。
【請求項7】
前記セグメンテーション手段は、任意に、前記固有空間上にサンプリング点を設定する請求項4又は請求項5記載の画像処理装置。
【請求項8】
前記目的関数は、前記特定の物体の形状の尤もらしさとして、前記形状パラメータが表す前記特定の物体の形状に対する画素における形状ラベルの値に関して単調に変化する単調関数を含む請求項1~請求項7の何れか1項記載の画像処理装置。
【請求項9】
前記探索アルゴリズムは、Branch and bound法及びグラフカット法である請求項3~請求項8の何れか1項記載の画像処理装置。
【請求項10】
受付手段、及びセグメンテーション手段を含み、特定の物体である被写体を表す入力画像から、前記被写体の領域を抽出する画像処理装置における画像処理方法であって、
前記受付手段が、前記入力画像を受け付けるステップと、
前記セグメンテーション手段が、(A)前記被写体を表す学習用の複数の画像であって、前記被写体の領域が予め求められた前記複数の画像に基づいて予め計算された固有ベクトルを基底とする固有空間であって、かつ、(B)固有空間上の点が、前記特定の物体の形状の統計的変動を表す統計的形状モデルの形状パラメータを示す固有空間において、前記入力画像に基づいて、前記固有空間上の点が示す前記形状パラメータが表す前記特定の物体の形状の尤もらしさと前記入力画像中の隣接画素間の画素値の差とに応じた値を表す予め定められた目的関数を最適化するように、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出するステップと、
を含む画像処理方法。
【請求項11】
前記セグメンテーション手段が、前記被写体の領域を、前記入力画像から抽出するステップは、
前記固有空間において、
前記入力画像に基づいて、前記目的関数を最適化するように、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定すると同時に、
前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出する請求項10記載の画像処理方法。
【請求項12】
前記セグメンテーション手段が、前記被写体の領域を、前記入力画像から抽出するステップは、
最適な形状パラメータを示す点が表す前記特定の物体の形状を含む形状集合を表す前記固有空間上の凸多胞体を繰り返し分割して最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索する探索アルゴリズムに従って、前記目的関数を最適化するように、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出する請求項11記載の画像処理方法。
【請求項13】
前記セグメンテーション手段が、前記被写体の領域を、前記入力画像から抽出するステップは、前記探索アルゴリズムにおいて、前記凸多胞体に含まれる形状集合に対する前記目的関数の下界を計算することにより、最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索して前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出する請求項12記載の画像処理方法。
【請求項14】
前記セグメンテーション手段が、前記被写体の領域を、前記入力画像から抽出するステップは、分割後の2つの凸多胞体の体積が対応するように、前記固有空間上の凸多胞体を分割し、かつ、最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索することを繰り返して、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出する請求項13記載の画像処理方法。
【請求項15】
前記セグメンテーション手段が、前記被写体の領域を、前記入力画像から抽出するステップは、前記固有空間上にサンプリング点を設定し、前記固有空間上に設定されたサンプリング点から決定される超平面を用いて前記凸多胞体を分割して最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索することを繰り返して、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、
前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出する請求項13又は請求項14記載の画像処理方法。
【請求項16】
前記セグメンテーション手段が、前記被写体の領域を、前記入力画像から抽出するステップは、任意に、前記固有空間上にサンプリング点を設定する請求項13又は請求項14記載の画像処理方法。
【請求項17】
前記目的関数は、前記特定の物体の形状の尤もらしさとして、前記形状パラメータが表す前記特定の物体の形状に対する画素における形状ラベルの値に関して単調に変化する単調関数を含む請求項11~請求項16の何れか1項記載の画像処理方法。
【請求項18】
前記探索アルゴリズムは、Branch and bound法及びグラフカット法である請求項12~請求項17の何れか1項記載の画像処理方法。
【請求項19】
コンピュータを、請求項1~請求項9の何れか1項記載の画像処理装置の各手段として機能させるためのプログラム。
発明の詳細な説明 【技術分野】
【0001】
本発明は、画像処理装置、方法、及びプログラムに係り、特に、被写体の領域を抽出する画像処理装置、方法、及びプログラムに関するものである。
【背景技術】
【0002】
デジタル画像の中からコンピュータを用いて図形を認識する処理はセグメンテーションと呼ばれるが、対象図形のSNが低く、かつ、形状が統計的に変動する場合、正確なセグメンテーションは非常に難しくなる。
【0003】
グラフカット(例えば、Boykov, Y., Veksler, O., Zabih, R.,“Fast approximate energy minimization via graph cuts.”, Pattern Analysis and Machine Intelligence, IEEE Transactions on 23 (11), 2001, p.1222‐1239.を参照。)などの最適化理論に基づくセグメンテーションアルゴリズムは、目的関数を真に最適化できるという意味で従来のアルゴリズムより優れており、SNが低い場合にもしばしば非常に高い性能を示した。しかし、それでも形状が変動する場合には正確に認識できないことがしばしばあった。
【0004】
そこで、最近は、最適化理論に基づくアルゴリズムにおいて図形の形状情報を利用する方法が主流になってきたが、それらは、大まかに次の3通りに分けられる。
【0005】
1つ目の方法としては、単一の形状テンプレートを利用する方法である。例えば、一般的形状テンプレートとして、楕円形状のテンプレート(例えば、Slabaugh, G., Unal, G., Sept, “Graph cuts segmentation using an elliptical shape prior.”, In: IEEE International Conference on Image Processing, 2005, Vol. 2. p.II‐1222‐5.を参照。)を用いる手法や、塊状図形のテンプレート(例えば、Funka-Lea, G., Boykov, Y., Florin, C., Jolly, M.‐P., Moreau-Gobard, R., Ramaraj, R., Rinck, D., 2006. “Automatic heart isolation for CT coronary visualization using graph-cuts.”, In: Biomedical Imaging: Nano to Macro, 3rd IEEE International Symposium on. IEEE, 2006, p.614‐617.を参照。)を利用する方法が知られている。
また、特定の形状テンプレートとして、ユーザ定義の任意形状(例えば、Freedman, D., Zhang, T., “Interactive graph cut based segmentation with shape priors.”, In: IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Vol. 1. IEEE, 2005, p.755‐762.を参照。)を用いる手法や、統計モデルから選択した形状(例えば、Grosgeorge, D., Petitjean, C., Dacher, J.-N., Ruan, S., “Graph cut segmentation with a statistical shape model in cardiac mri.”, Computer Vision and Image Understanding 117 (9), 2013, p.1027‐1035.,Akinobu Shimizu, Keita Nakagomi, Takuya Narihira, Hidefumi Kobatake, Shigeru Nawano, Kenji Shinozaki, Koich Ishizu, and Kaori Togashi, “Automated Segmentation of 3D CT Images based on Statistical Atlas and Graph Cuts”, Proc. of MICCAI workshop MCV, 2010, p.129-138.,Malcolm, J., Rathi, Y., Tannenbaum, A., “Graph cut segmentation with nonlinear shape priors.”, In: Image Processing, 2007. ICIP 2007. IEEE International Conference on. Vol. 4. IEEE, 2007, p.IV‐365.を参照。)を用いる手法が知られている。
【0006】
2つ目の方法としては、複数(数個)の形状テンプレートや形状の確率的表現を利用する方法である。例えば、統計モデルから選択した複数形状(例えば、Nakagomi, K., Shimizu, A., Kobatake, H., Yakami, M., Fujimoto, K., Togashi, K., “Multi-shape graph cuts with neighbor prior constraints and its application to lung segmentation from a chest CT volume.”, Medical image analysis 17 (1), 2013, p.62‐77.を参照。)を用いる手法や、形状の確率的表現(例えば、Linguraru, M. G., Pura, J. A., Pamulapati, V., Summers, R. M., “Statistical 4d graphs for multi-organ abdominal segmentation from multiphase CT.”, Medical image analysis 16 (4), 2012, p.904‐914.を参照。)を用いる手法が知られている。
【0007】
3つ目の方法としては、大量(数十以上)の形状テンプレートを利用する方法である。例えば、特定の問題用(例えば、Kohli, P., Rihan, J., Bray, M., Torr, P. H., “Simultaneous segmentation and pose estimation of humans using dynamic graph cuts.”, International Journal of Computer Vision 79 (3), 2008, p.285‐298.を参照。)のテンプレートを用いる手法や、任意の大量の図形形状(例えば、米国特許第8249349号明細書を参照。)を用いる手法が知られている。

【発明の概要】
【発明が解決しようとする課題】
【0008】
上述した従来の技術は、1つ目の手法から3つ目の手法の順に統計的に変動をする図形をうまく扱えるようになる。最も先端的な方法は、上記3つ目の手法の、任意の大量の図形形状(上記米国特許第8249349号明細書)である。しかし、この方法も含めてすべての方法は事前に形状を選択するなどして、あらかじめ形状集合を準備する必要がある。そのため、事前に選択した形状集合の中にそもそも最適な形状が含まれていなければ、セグメンテーションの最適性は保証されず、性能は低下する。
【0009】
一方、デジタル画像は有限個の画素の集合であることから、可能な全ての形状パターンを準備する方法も原理的には考えられ、その場合には最適性は保証される。しかし、そのために用意しなければならない形状数は膨大になり、前処理も含めた処理全体の必要なメモリと計算時間の観点から現実的ではない。例えば、従来最も優れていた上記米国特許第8249349号明細書のアルゴリズムでも、2次元画像で10個程度の数の形状テンプレートを扱うのが限界であった。しかし、統計モデルが生成しうる可能な形状パターン数は10個以上ある場合も多く、従来のアルゴリズムでは計算コストの観点から非常に困難である。さらに医用画像では3次元画像が主流であるが、その場合には必要なメモリは二桁以上増加し、ほぼ従来のアルゴリズムによって現実的な時間やメモリサイズで処理を実行することは不可能であった。
【0010】
本発明の一実施形態は、上記事情に鑑みてなされたものである。
【課題を解決するための手段】
【0011】
上記目的を達成するために、第1の態様の画像処理装置は、特定の物体である被写体を表す入力画像から、前記被写体の領域を抽出する画像処理装置であって、前記入力画像を受け付ける受付手段と、(A)前記被写体を表す学習用の複数の画像であって、前記被写体の領域が予め求められた前記複数の画像に基づいて予め計算された固有ベクトルを基底とする固有空間であって、かつ、(B)固有空間上の点が、前記特定の物体の形状の統計的変動を表す統計的形状モデルの形状パラメータを示す固有空間において、前記入力画像に基づいて、前記固有空間上の点が示す前記形状パラメータが表す前記特定の物体の形状の尤もらしさと前記入力画像中の隣接画素間の画素値の差とに応じた値を表す予め定められた目的関数を最適化するように、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出するセグメンテーション手段と、を備えている。
【0012】
また、第2の態様の画像処理方法は、受付手段、及びセグメンテーション手段を含み、特定の物体である被写体を表す入力画像から、前記被写体の領域を抽出する画像処理装置における画像処理方法であって、前記受付手段が、前記入力画像を受け付けるステップと、前記セグメンテーション手段が、(A)前記被写体を表す学習用の複数の画像であって、前記被写体の領域が予め求められた前記複数の画像に基づいて予め計算された固有ベクトルを基底とする固有空間であって、かつ、(B)固有空間上の点が、前記特定の物体の形状の統計的変動を表す統計的形状モデルの形状パラメータを示す固有空間において、前記入力画像に基づいて、前記固有空間上の点が示す前記形状パラメータが表す前記特定の物体の形状の尤もらしさと前記入力画像中の隣接画素間の画素値の差とに応じた値を表す予め定められた目的関数を最適化するように、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出するステップと、を含む。
【0013】
また、第3の態様の前記セグメンテーション手段は、前記固有空間において、前記入力画像に基づいて、前記目的関数を最適化するように、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定すると同時に、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出するようにすることができる。
【0014】
また、第4の態様の前記セグメンテーション手段は、最適な形状パラメータを示す点が表す前記特定の物体の形状を含む形状集合を表す前記固有空間上の凸多胞体を繰り返し分割して最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索する探索アルゴリズムに従って、前記目的関数を最適化するように、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出するようにすることができる。
【0015】
また、第5の態様の前記セグメンテーション手段は、前記探索アルゴリズムにおいて、前記凸多胞体の各頂点を調べることで、前記凸多胞体に含まれる形状集合に対する前記目的関数の下界を計算することにより、最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索して前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出するようにすることができる。
【0016】
また、第6の態様の前記セグメンテーション手段は、分割後の2つの凸多胞体の体積が対応するように、前記固有空間上の凸多胞体を分割し、かつ、最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索することを繰り返して、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出するようにすることができる。
【0017】
また、第7の態様の前記セグメンテーション手段は、前記固有空間上にサンプリング点を設定し、前記固有空間上に設定されたサンプリング点から決定される超平面を用いて前記凸多胞体を分割して最適な形状パラメータを示す点を含む前記固有空間上の凸多胞体を探索することを繰り返して、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出するようにすることができる。
【0018】
また、第8の態様の前記セグメンテーション手段は、任意に、前記固有空間上にサンプリング点を設定するようにすることができる。
【0019】
また、第9の態様の前記目的関数は、前記特定の物体の形状の尤もらしさとして、前記形状パラメータが表す前記特定の物体の形状に対する画素における形状ラベルの値に関して単調に変化する単調関数を含むようにすることができる。
【0020】
また、第10の態様の前記探索アルゴリズムは、Branch and bound法及びグラフカット法を用いることができる。
【0021】
また、第11の態様のプログラムは、コンピュータを、上記の画像処理装置の各手段として機能させるためのプログラムである。
【発明の効果】
【0022】
実施の形態に係る画像処理装置、方法、及びプログラムによれば、(A)被写体の領域が予め求められた複数の画像に基づいて予め計算された固有ベクトルを基底とする固有空間であって、(B)固有空間上の点が、特定の物体の形状を表す統計的形状モデルの形状パラメータを示す固有空間において、固有空間上の点が示す形状パラメータが表す特定の物体の形状の尤もらしさと入力画像中の隣接画素間の画素値の差とに応じた値を表す予め定められた目的関数を最適化するように、入力画像が表す被写体の形状を表す形状パラメータを推定する。推定された形状パラメータが表す特定の物体の形状を事前知識として、被写体の領域を、入力画像から抽出する。これにより、計算量の増大を抑制して、被写体の領域を精度よく抽出することができる、という効果が得られる。
【図面の簡単な説明】
【0023】
【図1】実施の形態に係る統計的形状モデル生成装置の構成を示すブロック図である。
【図2】実施の形態に係る統計的形状モデル生成装置及び画像処理装置のコンピュータ構成例を示すブロック図である。
【図3】肝臓の3次元画像例を示す説明図である。
【図4】肝臓の形状の画像を対象として主成分分析を行った場合の、第1主成分に沿った肝臓の形状の統計的ばらつきを示す図である。
【図5】統計的形状モデルの例を示す説明図である。
【図6】固有空間と形状空間との関係を示す説明図である。
【図7】実施の形態に係る画像処理装置の構成を示すブロック図である。
【図8】固有空間の分割法を説明するための説明図である。
【図9】実施の形態における学習処理ルーチンを示すフローチャートである。
【図10】実施の形態における画像処理ルーチンを示すフローチャートである。
【図11】実施の形態におけるセグメンテーション処理ルーチンを示すフローチャートである。
【図12】格子状にサンプリング点を設定する場合の近似解法を説明するための図である。
【図13】第2の実施の形態における140例の膵臓の認識結果を示す図である。
【図14】第2の実施の形態における計算時間の比較結果を示す図である。
【図15】ランダムにサンプリング点を設定する場合の近似解法を説明するための図である。
【図16】目的関数の下界が単調性を満たす関係を示す図である。
【図17】目的関数の下界が単調性を満たすh(y)とh(y)との一例を示す図である。
【図18】統計的形状モデルとしてLogOddsを用いる場合の例を示す図である。
【図19】従来技術を示す説明図である。
【発明を実施するための形態】
【0024】
以下、図面を参照して、実施の形態を詳細に説明する。なお、実施の形態では、特定の物体の統計的変動を表す図形の形状の統計モデル(以下、統計的形状モデルと称する。)を生成する統計的形状モデル生成装置10と、特定の物体である被写体を表す入力画像から、被写体の領域を抽出する画像処理装置100とについて説明する。また、実施の形態では、3次元腹部CT画像を入力画像とし、当該入力画像から、膵臓の領域を抽出する場合を例に説明する。

【0025】
<概要>
本実施の形態では、統計的形状モデルを用いるセグメンテーションアルゴリズムを提案する。このアルゴリズムの特徴は以下の5点である。

【0026】
1.セグメンテーションの過程で、統計的形状モデルを用いて生成可能な全ての3次元形状(10個以上)を考慮可能であり、セグメンテーションの目的関数の観点から最適な形状を選択可能である。

【0027】
2.大量の形状テンプレートを用いる従来手法では必須で、特に計算コストが高い前処理(形状の生成・選択とクラスタリング)が不要である。

【0028】
3.Branch and bound法とグラフカット法との組み合わせに限らず、最適化理論に基づくアルゴリズムであれば適用可能である。

【0029】
4.異なる統計的形状モデルや固有空間の分割法に対しても適用可能である。

【0030】
5.3次元腹部CT画像内の膵臓のセグメンテーションに関して世界最高の精度が得られる。

【0031】
図19に従来技術の処理の概要を示す。図19は、上記米国特許第8249349号明細書の手法を、統計的形状モデルを利用した膵臓のセグメンテーションに適用させた場合の処理である。図19に示すように、従来手法では、形状テンプレート集合T(⊂S)の生成処理、及び集合Tに対するクラスタリング処理が必要であるが、本実施の形態では、上述したように、大量の形状テンプレートを用いる従来手法で必要であった前処理が不要である。なお、従来手法(上記米国特許第8249349号明細書)で扱われていた画像は2次元であり、形状の最大数は10であるのに対し、本実施の形態で用いるアルゴリズムは3次元形状(10個以上)を考慮可能である。

【0032】
また、実施の形態では、統計的形状モデルとしてレベルセット分布モデル(Level set distribution model; LSDM。符号付距離モデル(Signed distance model)と呼ばれることもある)を利用する。以下では、まず、統計的形状モデル生成装置10において統計的形状モデル(LSDM)の固有ベクトルと固有値とを算出した後で、画像処理装置100において統計的形状モデルの形状パラメータを推定し、推定された形状パラメータを用いた統計的形状モデルに基づいた最適セグメンテーションアルゴリズムを示す。

【0033】
[第1の実施の形態]
<第1の実施の形態に係る統計的形状モデル生成装置10の構成>
第1の実施の形態に係る統計的形状モデル生成装置10は、図1に示される機能ブロックで表すことができる。また、これらの機能ブロックは、図2に示されるコンピュータのハードウェア構成により実現することができる。図2を参照してコンピュータの構成を説明する。

【0034】
図2に示す第1の実施の形態に係る統計的形状モデル生成装置10は、プログラムに基づき統計的形状モデル生成装置10の本実施の形態に係る処理を行うCPU(Central Processing Unit;中央処理装置)21と、CPU21による各種プログラムの実行時のワークエリア等として用いられるRAM(Random Access Memory)22と、各種制御プログラムや各種パラメータ等が予め記憶された記録媒体であるROM(Read Only Memory)23と、各種情報を記憶するために用いられるハードディスク24(図中「HDD」と記載)と、キーボードやマウス等からなる入力装置25と、ディスプレイ等からなる表示装置26と、LAN(Local Area Network)等を用いて通信を行う通信装置27と、外部に接続された画像情報提供装置30との間の各種情報の授受を司る入出力インタフェース部(図中、「外部IF」と記載)28と、を備えており、これらがシステムバスBUS29により相互に接続されて構成されている。

【0035】
CPU21は、RAM22、ROM23、及びハードディスク24に対するアクセス、入力装置25を介した各種情報の取得、表示装置26に対する各種情報の表示、通信装置27を用いた各種情報の通信処理、及び入出力インタフェース部28に接続された画像情報提供装置30を含む外部装置からの情報の入力等を、各々行うことができる。

【0036】
CPU21が、ハードディスク24に記憶された本実施形態に係る処理を制御するプログラムを、RAM22に読み込み実行することにより、図1に示す本実施の形態に係る統計的形状モデル生成装置10における、図1に示す各処理部の機能が実行される。

【0037】
このようなコンピュータ構成により、図1に示す本実施の形態に係る統計的形状モデル生成装置10が構成されている。なお、図1は機能ブロックとなる構成を表し、一方、図2はデバイス等の接続状態を表すものである。前記したように、機能ブロックとデバイス等とは有機的、かつ相互に関連して統計的形状モデル生成装置10を構成するため、図2と共に図1についても詳細に説明する。

【0038】
統計的形状モデル生成装置10は、図2に示したコンピュータのハードウェア及び制御プログラムを含むソフトウェアを利用して構成される機能として、学習用受付部12と、学習部14と、統計的形状モデルデータベース16とを備えている。

【0039】
学習用受付部12は、特定の物体である被写体を表す複数の画像であって、被写体の領域が予め求められた複数の画像を学習データとして受け付ける。本実施の形態で被写体となる特定の物体は膵臓であるため、学習用受付部12は、膵臓の領域が予め求められた複数の画像を受け付ける。

【0040】
学習部14は、学習用受付部12によって受け付けた膵臓の領域が予め求められた複数の画像に基づいて、固有ベクトルと固有値とを算出する。例えば、学習部14は、膵臓の領域が予め求められた複数の画像に基づいて、主成分分析によって、固有ベクトルと固有値とを算出する。固有ベクトルが算出されることにより、固有ベクトルを基底とする固有空間が生成される。また、学習部14は、算出された固有ベクトルと固有値とを統計的形状モデルデータベース16に格納する。

【0041】
ここで、図3に、肝臓の形状の一例を示す(なお、図3に示す例は肝臓であるが、問題の本質は膵臓の場合と同様である。)。図3に示すように、肝臓の形状は様々である。そのため、実施の形態では、特定の物体の形状を表す統計的形状モデルを用いて、形状のばらつきを少数の形状パラメータにより表現する。また、図4に、主成分分析によって得られた第1主成分に対応する肝臓の形状の固有ベクトルを示す。図4に示すように、固有ベクトルに対応する形状パラメータαの値に応じて、表わされる肝臓の形状が変化することがわかる。なお、σは、統計的ばらつきを表すパラメータである。

【0042】
図5に、統計的形状モデルの一例を示す。図5に示すように、統計的形状モデルとしては様々なものが存在するが、本実施の形態では、統計的形状モデルとして、LSDM(例えば、参考文献(Cremers, D., Rousson, M., Deriche, R., “A review of statistical approaches to level set segmentation: integrating color, texture, motion and shape.”, International journal of computer vision 72 (2), 2007, p.195‐215.)を参照。)を用いる。

【0043】
LSDMは、代表的な統計的形状モデルであり、学習データの図形ラベルに対する符号付の距離画像を作成し、符号付の距離画像に対して線形の統計解析(例えば、PCAやICAなど)を行うことで任意形状を以下の式(1)で表現する方法である。

【0044】
【数1】
JP2016027840A1_000003t.gif

【0045】
ここで、上記式(1)は画素p(∈P)に関する式である。なお、Pは画素集合を表し、φ(α)は、画素pに関する任意形状に対応するレベルセット関数を表し、μは、画素pに関する平均のレベルセット関数を表し、λはi番目の固有値を表し、{u,...,u}はd個の固有ベクトルにおける画素pに対する成分を表し、α=[α,... ,αは、固有空間上のRα={r∈R|||r||≦w}の領域(通常は±3σの範囲)で定義される形状パラメータを表す。また、wは、正の定数を示す。デジタル画像における形状は、パラメータαを以下の式(2)に示す関数gで写像することによって得られる。なお、実施の形態では、固有空間の次元d=2である場合を例に説明する。

【0046】
【数2】
JP2016027840A1_000004t.gif

【0047】
ここで、L={0,1}はラベル集合で、0は背景であることを示し、1は図形であることを示す。yは画像における1つの形状を表す。また、LSDMの場合、以下の式(3)に示すように、パラメータαを形状に写像する関数gは、Heaviside function H(・)である。

【0048】
【数3】
JP2016027840A1_000005t.gif

【0049】
図6に、d=2である場合のLSDMの固有空間(固有ベクトルが張る空間)とデジタル画像における形状集合Sとの関係を示す。関数gは固有空間からデジタル画像の形状集合S上への写像である。また、形状集合S⊂L|P|はgの像である。図6に示すように、固有空間では形状ラベルの集合は多角形を構成する。各形状y∈Sは固有空間内に原像を持つ。

【0050】
統計的形状モデルデータベース16には、学習部14によって算出された固有ベクトルと固有値とが格納される。

【0051】
<画像処理装置100の構成>
本実施の形態に係る画像処理装置100は、図7に示される機能ブロックで表すことができる。また、これらの機能ブロックは、上記図2に示されるコンピュータのハードウェア構成により実現することができる。

【0052】
画像処理装置100は、上記図2に示したコンピュータのハードウェア及び制御プログラムを含むソフトウェアを利用して構成される機能として、図7に示すように、受付部102と、演算部104と、出力部120とを備えている。

【0053】
受付部102は、入力画像として、3時相の3次元腹部CT画像を受け付ける。なお、受付部102は、早期相画像、門脈相画像、及び晩期相画像を、3時相の3次元腹部CT画像として受け付ける。

【0054】
演算部104は、受付部102によって受け付けた3次元腹部CT画像から、膵臓の領域を抽出する。演算部104は、統計的形状モデルデータベース106と、画像処理部108とを備えている。

【0055】
統計的形状モデルデータベース106には、統計的形状モデルデータベース16と同じ固有ベクトルと固有値とが格納されている。

【0056】
画像処理部108は、画像間位置合わせ部110と、空間的標準化部112と、セグメンテーション部114とを備えている。

【0057】
画像間位置合わせ部110は、受付部102によって受け付けた、早期相画像と、門脈相画像と、晩期相画像との間の位置合わせを行い(例えば、Shimizu, A., Kimoto, T., Kobatake, H., Nawano, S., Shinozaki, K., “Automated pancreas segmentation from three-dimensional contrast-enhanced computed tomography.”, International journal of computer assisted radiology and surgery 5 (1), 2010, p.85‐98.を参照。)、位置合わせ画像を生成する。

【0058】
空間的標準化部112は、画像間位置合わせ部110によって生成された位置合わせ画像に基づいて、例えば予め定められた非線形関数を用いて、位置合わせ画像を標準となる画像に合わせるように空間的標準化を行い(例えば、Shimizu, A., Kimoto, T., Kobatake, H., Nawano, S., Shinozaki, K., “Automated pancreas segmentation from three-dimensional contrast-enhanced computed tomography.”, International journal of computer assisted radiology and surgery 5 (1), 2010, p.85‐98.を参照。)、空間的標準化画像を生成する。

【0059】
セグメンテーション部114は、統計的形状モデルデータベース106に格納された予め計算された固有ベクトルを基底とする固有空間において、受付部102によって受け付けた入力画像に基づいて、目的関数を最適化するように、空間的標準化部112によって生成された空間的標準化画像が表す被写体の形状を表す形状パラメータを推定し、空間的標準化画像が表す被写体の領域を求める。具体的には、セグメンテーション部114は、形状パラメータを推定すると同時に、被写体の領域を抽出する。

【0060】
なお、上記図6に示すように、固有ベクトルを基底とする固有空間上の点は、統計的形状モデルの形状パラメータを示す。また、目的関数は、固有空間上の点が示す形状パラメータが表す特定の物体の形状の尤もらしさと入力画像中の隣接画素間の画素値の差とに応じた値を表すように、予め定められている。

【0061】
一般に、上述した形状集合L|P|のサイズ(=デジタル画像上の形状数)は膨大であり、従来のアルゴリズムでは扱うことができなかった。本実施の形態におけるアルゴリズムでは、目的関数の最適化を固有空間において実行する点が従来とは異なり、これによって膨大な数の形状を効率的に扱うことが可能である。

【0062】
セグメンテーション部114は、固有空間上の凸多角形を繰り返し分割して最適な形状パラメータを示す点を含む凸多角形を探索する探索アルゴリズムに従って、以下の式(4)の目的関数を最適化するように、入力画像が表す被写体の形状を表す形状パラメータを推定する。なお、探索される固有空間上の凸多角形は、最適な形状パラメータを示す点が表す特定の物体の形状を表している。
また、セグメンテーション部114は、探索アルゴリズムにおいて、凸多角形の各頂点を調べることで、その凸多角形に含まれる形状集合に対する目的関数の下界を計算することにより、最適な形状パラメータを示す点を含む固有空間上の凸多角形を探索する。
なお、本実施の形態では、固有空間が2次元である場合を例に説明するため、セグメンテーション部114は固有空間上の凸多角形を探索するが、固有空間が3次元以上の場合には、セグメンテーション部114は固有空間上の凸多胞体を探索する。また、探索時に固有空間を分割する直線は、固有空間が3次元以上の場合には超平面となる。

【0063】
本実施の形態における重要なアイディアは、上記図6に示したように、デジタル画像上の形状集合が、固有空間上の凸多角形に対応する事実を利用することにある。本発明の実施の形態では、当該事実を利用することにより、多数の形状を一度に効率的に扱うことが可能である。

【0064】
セグメンテーション部114は、デジタル画像上の複数の形状集合が、固有空間上の凸多角形に対応する事実に基づいて、予め設定した目的関数を最小化する形状に対応する、固有空間上の凸多角形を探索する。そして、セグメンテーション部114は、探索された固有空間上の凸多角形に基づいて、目的関数を最小化する形状を形状集合Sから見つけ出し、見つけ出された形状を事前知識として当該形状を被写体の領域として抽出する。本実施の形態では、グラフカット(上記Boykov, Y., Veksler, O., Zabih, R.,“Fast approximate energy minimization via graph cuts.”, Pattern Analysis and Machine Intelligence, IEEE Transactions on 23 (11), 2001, p.1222‐1239.を参照。)で良く用いられるエネルギー関数を、予め設定した目的関数として用いる場合を例に説明する。本実施の形態で用いる目的関数を、以下の式(4)に示す。

【0065】
【数4】
JP2016027840A1_000006t.gif

【0066】
ここで、Nは、隣接する画素のペアの集合である。xは画素pのラベルを表し、背景であれば0の値をとり、図形であれば1の値をとる。また、Iは、画素pの画素値を表し、yは、画素pにおける推定形状の値を表し、推定形状が背景であれば0の値をとり、図形であれば1の値をとる。

【0067】
また、上記式(4)におけるF(I,y)と、B(I,y)とは、以下の式(5)、(6)で定義される。ここで、以下の式(5)は、画素pが図形として割り当てられるコストを表し、以下の式(6)は、画素pが背景として割り当てられるコストを表す。

【0068】
【数5】
JP2016027840A1_000007t.gif

【0069】
は画素pの画素値を表し、Iは画素qの画素値を表す。また、λ、λは予め定められた正の定数である。また、上記式(4)におけるPpq(I,I)は、以下の式(7)で定義され、互いに隣接する画素pと画素qの間の画素値の差を評価している。

【0070】
【数6】
JP2016027840A1_000008t.gif

【0071】
本実施の形態では、上記式(4)に示す目的関数の最適化手法として、一般的な最適化アルゴリズムの一つであるbranch and bound探索アルゴリズムを用いる。本実施の形態では、形状集合に対する目的関数の下界を効率的に計算する方法を提案する。形状集合に対する目的関数の下界を効率的に計算するための重要なアイディアは、凸多角形の頂点のみを調べれば、その凸多角形に含まれるすべての形状集合に対する目的関数の下界が分かる点である。

【0072】
上記式(4)を、形状yの代わりに、関数g(α)を用いて表すと、以下の式(8)で表すことができる。

【0073】
【数7】
JP2016027840A1_000009t.gif

【0074】
[branching処理]
branch and bound探索アルゴリズムのbranching処理では、親ノードH(∈Rα)が与えられると、親ノードHは子ノードH及びHへ分解され、上記式(8)に示した目的関数は、以下の式(9)で表わされる。

【0075】
【数8】
JP2016027840A1_000010t.gif

【0076】
ここで、branching処理における分割処理は、様々な分割方法が可能である。本実施の形態では、上記式(1)に示したφ(α)の符号によって、以下の式(10)及び(11)に示すように、親ノードHは子ノードH及びHへ分割される。なお、画素kは、集合Qからサンプリングによって選択された画素を表す。

【0077】
【数9】
JP2016027840A1_000011t.gif

【0078】
ここで、集合Qは、以下の式(12)に示すように、固有空間内のノードHを切る直線φ(α)=0であるような画素kの集合である。なお、固有空間が3次元以上である場合、集合Qは、固有空間内のノードHを切る超平面φ(α)=0であるような画素kの集合となる。

【0079】
【数10】
JP2016027840A1_000012t.gif

【0080】
また、Vは、ノードHが表す凸多角形の頂点集合を表す。Vは、Hのエッジに対応するφ(α)の連立方程式を解析的に解くことで得ることができる。

【0081】
[bounding処理]
branch and bound探索アルゴリズムのbranching処理においては、下界L(H)(i={1,2})は、以下の式(13)~(14)で与えられる。

【0082】
【数11】
JP2016027840A1_000013t.gif

【0083】
上記式(13)は、最小化問題のためのJensenの不等式である。上記式(14)における「max」と「min」との変化は、上記式(5)のマイナス符号に基づいている。また、α∈Hにおけるφ(α)の最大値と最小値とは、線形計画法の基礎理論に基づき、Hの頂点集合Vから取得され、以下の式(15)及び(16)で表すことができる。

【0084】
【数12】
JP2016027840A1_000014t.gif

【0085】
また、セグメンテーション部114は、分割後の2つの凸多角形の面積が対応するように(面積がほぼ等しくなるように)、固有空間上の凸多角形を分割することを繰り返すことにより、最適な形状パラメータを示す点を含む固有空間上の凸多角形を探索する。

【0086】
セグメンテーション部114によるセグメンテーションアルゴリズムの疑似コードを表1に示す。まず、対象画像Iが与えられると、セグメンテーション部114は、固有空間全体Rα={r∈R|||r||≦w}を親ノードHとして処理を開始する。すなわち、セグメンテーション部114は、Rαをルートとして処理を開始し、デジタル画像から関数Select(Q)を用いて選択した画素kによって親ノードHを分割する。分割後の子ノードHとHとは、Queueに入れられる。ただし、Queueに入れられたノードは、下界の値が昇順になるように並べられているとする。次に、セグメンテーション部114は、Queueから、最も低い下界のノードを選択して親ノードHとし、上記の処理をQが空集合になるまで反復する。反復終了時の結果から、最適な形状パラメータが得られ、最適な形状が求まる。そして、セグメンテーション部114は、最適な形状を用いて、最適なセグメンテーション結果を得る。

【0087】
【表1】
JP2016027840A1_000015t.gif

【0088】
なお、Hは、固有空間における凸多角形を表す。なお、Hの初期値はRαに対応する。また、H、Hは、Hを分割してできた二つの凸多角形を表す。Qは、Hを分割可能なデジタル画像上の画素集合である。Queueは、凸多角形ノードが入っているキュー(Branch and bound探索アルゴリズムで利用)を表す。Select(Q)は、分割のための画素kを選択する関数を表す。

【0089】
また、効率的にアルゴリズムを動作させるためのSelect(Q)に関する新しい提案も、本実施の形態における重要な特徴である。ここで、Select(Q)は画像から凸多角形を分割するための画素を選択する関数であるが、探索を効率的に行なうために、分割後の子ノードH,Hに対応する凸多角形の面積がほぼ等しくなるような分割法を本実施の形態で提案し、その効率を実際の画像で検証した。

【0090】
本実施の形態では、以下の式(17)に示す関数を、Select(Q)として用いる。

【0091】
【数13】
JP2016027840A1_000016t.gif

【0092】
ここで、

【0093】
【数14】
JP2016027840A1_000017t.gif

【0094】
である。

【0095】
図8に、上記式(17)の関数の作用の概念図を示す。図8に示すように、上記式(17)の作用は、直感的には、Hを面積が等しい2つの領域(H,H)に分割するような画素pを選択する。

【0096】
セグメンテーション部114は、探索された特定の物体の最適な形状パラメータが表す形状を事前知識として、目的関数を最適化するように、特定の物体の形状である被写体の領域を、空間的標準化部112によって生成された空間的標準化画像から抽出する。

【0097】
出力部120は、セグメンテーション部114によって抽出された被写体の領域を結果として出力する。

【0098】
<統計的形状モデル生成装置10の作用>
次に、統計的形状モデル生成装置10の作用を説明する。

【0099】
図9に示す統計的形状モデル生成装置10の処理は、上記図2におけるCPU21によりHDD24等に記憶されたプログラムに基づき行なわれるものである。統計的形状モデル生成装置10に、膵臓の領域が予め求められた複数の画像が入力されると、図9に示す学習処理ルーチンが実行される。

【0100】
まず、ステップS100において、学習用受付部12は、膵臓の領域が予め求められた複数の画像を受け付ける。

【0101】
次に、ステップS102において、学習部14は、上記ステップS100で受け付けた膵臓の領域が予め求められた複数の画像に基づいて、主成分分析によって、固有ベクトルと固有値とを算出する。

【0102】
ステップS104において、学習部14は、上記ステップS102で算出された固有ベクトルと固有値とを統計的形状モデルデータベース16に格納し、学習処理ルーチンを終了する。

【0103】
<画像処理装置100の作用>
次に、画像処理装置100の作用を説明する。

【0104】
図10に示す画像処理装置100の処理は、上記図2におけるCPU21によりHDD24等に記憶されたプログラムに基づき行なわれるものである。まず、統計的形状モデル生成装置10の統計的形状モデルデータベース16に格納されている固有ベクトルと固有値とが、画像処理装置100に入力されると、統計的形状モデルデータベース106に格納される。そして、画像処理装置100に、3時相の3次元腹部CT画像が入力画像として入力されると、上記図10に示す画像処理ルーチンが実行される。

【0105】
ステップS200において、受付部102は、入力画像として、3時相の3次元腹部CT画像を受け付ける。

【0106】
ステップS202において、画像間位置合わせ部110は、上記ステップS200で3時相の3次元腹部CT画像として受け付けた、早期相画像と、門脈相画像と、晩期相画像との間の位置合わせを行い、位置合わせ画像を生成する。

【0107】
ステップS204において、空間的標準化部112は、上記ステップS202で生成された位置合わせ画像に基づいて、例えば非線形関数を用いて、空間的標準化画像を生成する。

【0108】
ステップS206において、セグメンテーション部114は、統計的形状モデルデータベース106に格納された予め計算された固有ベクトルを基底とする固有空間において、上記ステップS200で受け付けた入力画像に基づいて、目的関数を最適化するように、上記ステップS204で生成された空間的標準化画像が表す被写体の形状を表す形状パラメータを推定し、空間的標準化画像が表す被写体の形状を求める。ステップS206は、図11に示すセグメンテーション処理ルーチンによって実現される。

【0109】
<セグメンテーション処理ルーチン>
ステップS300において、セグメンテーション部114は、固有空間全体Rαを親ノードHとして設定する。

【0110】
ステップS302において、セグメンテーション部114は、上記式(12)に示すように、親ノードHを分割可能な画素kの集合を、集合Qとして設定する。

【0111】
ステップS304において、セグメンテーション部114は、Queueを初期化する。

【0112】
ステップS306において、セグメンテーション部114は、上記ステップS300で設定した親ノードHと、上記式(14)に従って算出される当該親ノードHに対応する下界L(H;I)との組み合わせを、上記ステップS304で初期化されたQueueに格納する。

【0113】
ステップS308において、セグメンテーション部114は、上記ステップS302で設定された集合Q、又は後述するステップS316で更新された集合Qから、上記式(17)に従って、画素kを選択する。

【0114】
ステップS310において、セグメンテーション部114は、上記ステップS300で設定された親ノードH、又は後述するステップS314で更新された親ノードHを、上記ステップS308で選択された画素kを用いて上記式(10)及び(11)に従って分割し、子ノードH及びHを設定する。

【0115】
ステップS312において、セグメンテーション部114は、上記ステップS310で設定された子ノードH及び上記式(14)に従って算出される当該子ノードHに対応する下界L(H;I)との組み合わせを、Queueに格納する。また、上記ステップS310で設定された子ノードH及び上記式(14)に従って算出される当該子ノードHに対応する下界L(H;I)との組み合わせを、Queueに格納する。

【0116】
ステップS314において、セグメンテーション部114は、Queueに格納されたノードのうち、最も低い下界のノードを選択し、親ノードHを、選択されたノードに更新する。

【0117】
ステップS316において、セグメンテーション部114は、上記ステップS314で更新された親ノードHに基づいて、上記式(12)に示すように、集合Qを、親ノードHを分割可能な画素kの集合に更新する。なお、親ノードHを分割可能な画素kが存在しない場合には、集合Qが空集合に更新される。

【0118】
ステップS318において、セグメンテーション部114は、集合Qが空集合であるか否かを判定する。集合Qが空集合である場合には、ステップS308へ戻る。一方、集合Qが空集合でない場合には、ステップS320へ進む。

【0119】
ステップS320において、セグメンテーション部114は、上記ステップS314で最終的に更新された親ノードHに含まれる形状パラメータαを選択し、選択された形状パラメータαを関数gに代入し、特定の物体の形状yを決定する。

【0120】
ステップS322において、セグメンテーション部114は、上記ステップS320で決定された特定の物体の形状y及び目的関数に基づいて、目的関数を最適化するように、特定の物体の形状である被写体の領域xを、上記ステップS204で生成された空間的標準化画像から抽出する。

【0121】
ステップS324において、セグメンテーション部114は、上記ステップS322で抽出された被写体の領域xを結果として出力する。

【0122】
次に画像処理ルーチンに戻り、出力部120は、ステップS208において、上記ステップS206で抽出された被写体の領域xを結果として出力して、画像処理ルーチンを終了する。

【0123】
以上説明したように、本実施形態の画像処理装置100では、プログラムに基づくコンピュータ処理による機能として、(A)被写体の領域が予め求められた複数の画像に基づいて予め計算された固有ベクトルを基底とする固有空間であって、(B)固有空間上の点が、特定の物体の形状を表す統計的形状モデルの形状パラメータを示す固有空間において、固有空間上の点が示す形状パラメータが表す特定の物体の形状の尤もらしさと入力画像中の隣接画素間の画素値の差とに応じた値を表す予め定められた目的関数を最適化するように、入力画像が表す被写体の形状を表す形状パラメータを推定する。推定された形状パラメータが表す特定の物体の形状を事前知識として、被写体の領域を、入力画像から抽出する。これにより、計算量の増大を抑制して、被写体の領域を精度よく抽出することができる。

【0124】
また、本実施形態の画像処理装置100は、被写体の領域を高速に抽出することができる。

【0125】
また、提案するアルゴリズムと従来のアルゴリズムとの重要な違いの一つに、前処理が不要である点が挙げられる。従来は、あらかじめ形状を用意し、クラスタリングをする前処理が必要であり、前処理にかなりの時間とメモリを必要としていたが、本実施の形態における提案アルゴリズムでは前処理は一切不要である。これにより、コストの削減を実現した。

【0126】
また、本実施形態の形態では、140例の3次元腹部CT画像から膵臓をセグメンテーションする実験を行って処理精度や処理時間を検証し、以下の結果が得られた。

【0127】
1.3次元腹部CT画像からの膵臓のセグメンテーションに関して世界最高の精度が得られた。
2.提案アルゴリズムによる目的関数の最小化後の値は、従来手法と比較すると必ず小さくなることが理論的に保証されているが、そのことを実際の画像でも確認した。
3.従来手法(上記米国特許第8249349号明細書)で扱われていた形状の最大数が10であるのに比べ、提案法では10以上の数を扱うことができるようになった。また、画像も従来は2次元の画像に限られていたのに対して、画素数が二桁近く多い3次元画像も扱えるようになった。

【0128】
[第2の実施の形態]
次に、第2の実施の形態について説明する。なお、第2の実施の形態に係る統計的形状モデル生成装置及び画像処理装置の構成は、第1の実施の形態と同様の構成となるため、同一符号を付して説明を省略する。

【0129】
第2の実施の形態では、近似解法を用いて、被写体の形状を表す形状パラメータを推定する点が、第1の実施の形態と異なる。

【0130】
図12に、第2の実施の形態における近似解法を説明するための図を示す。図12に示すように、上記第1の実施の形態で用いた厳密解法とは異なり、第2の実施の形態では近似解法によって形状パラメータを推定する。第2の実施の形態では、図12に示すように、固有空間上に格子状にサンプリング点を設定し、サンプリング点から決定される直線を用いて凸多角形を分割することを繰り返す。なお、サンプリング点の設定は、予め定められたサンプリングパラメータkに応じて行われる。例えば、サンプリングパラメータk=4の場合、固有空間の次元ごとに2個のサンプリング点が設定される。

【0131】
第2の実施の形態に係る画像処理装置のセグメンテーション部114は、まず、固有空間上に格子状にサンプリング点を設定する。

【0132】
そして、セグメンテーション部114は、固有空間上に設定されたサンプリング点から、分割後の2つの凸多角形の面積が対応するように、直線を設定し、設定された直線で凸多角形を分割して、最適な形状パラメータを示す点を含む凸多角形を探索することを繰り返し、入力画像が表す被写体の形状を表す形状パラメータを推定する。
具体的には、上記表1に示したセグメンテーションアルゴリズムの疑似コードにおいて、関数Select(Q)を用いて選択した画素kによって親ノードHを分割する処理の代わりに、サンプリング点から決定される直線を用いて、親ノードHを分割する処理を行う。

【0133】
第2の実施の形態では、近似解法を用いることにより、固有空間の次元数を増加させることができ、形状パラメータを精度良く推定することができる。

【0134】
図13に、第2の実施の形態に係る画像処理装置による140例の膵臓の認識結果を示す。図13において、破線の円は上記第1の実施の形態に係る画像処理装置による膵臓の認識結果を表し、実線の円は第2の実施の形態に係る画像処理装置による膵臓の認識結果を表す。
図13に示すように、近似解法を用いることにより、固有空間の次元数を上げることで厳密解よりも高いセグメンテーション精度が達成されていることがわかる。

【0135】
また、図14に、第2の実施の形態の近似解法、上記米国特許第8249349号明細書に記載の手法を膵臓のセグメンテーションに適用させた場合の解法、及び上記第1の実施の形態の厳密解法の計算時間の比較結果を示す。図14は、固有空間の次元dを増加させた場合の計算時間の変化を示している。なお、厳密解法については計算時間の上限を100hとした。また、近似解法についてはサンプリングパラメータの値をk=4とした。また、CPUは、3.1GHzのIntel(R)Xeon(R)を用い、前処理に1CPU、最適化に2CPUを用いた。

【0136】
図14に示すように、近似解法を用いた場合には、上記第1の実施の形態の厳密解法を用いた場合及び上記米国特許第8249349号明細書に記載の手法を膵臓のセグメンテーションに適用させた場合に比べ、固有空間の次元数を上げた場合の計算コストを大幅に削減可能であることがわかる。

【0137】
なお、第2の実施の形態に係る画像処理装置の他の構成及び作用については、第1の実施の形態と同様であるため、説明を省略する。

【0138】
以上説明したように、第2の実施の形態に係る画像処理装置によれば、固有空間上の凸多角形上に、格子状にサンプリング点を設定し、固有空間上の凸多角形に設定された直線を用いて凸多角形を分割して探索することを繰り返して、入力画像が表す被写体の形状を表す形状パラメータを推定する。推定された形状パラメータが表す特定の物体の形状を事前知識として、被写体の領域を、入力画像から抽出する。これにより、計算量の増大を抑制して、被写体の領域を精度よく抽出することができる。

【0139】
また、近似解法により、次元数を増加させたより高度な統計的形状モデルが利用可能となる。

【0140】
[第3の実施の形態]
次に、第3の実施の形態について説明する。なお、第3の実施の形態に係る統計的形状モデル生成装置及び画像処理装置の構成は、第1の実施の形態と同様の構成となるため、同一符号を付して説明を省略する。

【0141】
第3の実施の形態では、ランダムにサンプリング点を設定する近似解法を用いて、被写体の形状を表す形状パラメータを推定する点が、第2の実施の形態と異なる。

【0142】
図15に、第3の実施の形態における近似解法を説明するための図を示す。図15に示すように、第3の実施の形態では、固有空間上にランダムにサンプリング点を設定し、サンプリング点から決定される直線を用いて凸多角形を分割することを繰り返す。

【0143】
第3の実施の形態に係る画像処理装置のセグメンテーション部114は、まず、固有空間上にランダムにサンプリング点を設定する。

【0144】
そして、セグメンテーション部114は、固有空間上に設定されたサンプリング点から、分割後の2つの凸多角形の面積が対応するように、直線を設定し、設定された直線で凸多角形を分割して、最適な形状パラメータを示す点を含む凸多角形を探索することを繰り返し、入力画像が表す被写体の形状を表す形状パラメータを推定する。

【0145】
具体的には、上記表1に示したセグメンテーションアルゴリズムの疑似コードにおいて、関数Select(Q)を用いて選択した画素kによって親ノードHを分割する処理の代わりに、サンプリング点から決定される直線を用いて親ノードHを分割する処理を行う。

【0146】
なお、第3の実施の形態に係る画像処理装置の他の構成及び作用については、第1又は第2の実施の形態と同様であるため、説明を省略する。

【0147】
以上説明したように、第3の実施の形態に係る画像処理装置によれば、固有空間上の凸多角形上に、ランダムにサンプリング点を設定し、固有空間上の凸多角形に設定された直線を用いて凸多角形を分割して探索することを繰り返して、入力画像が表す被写体の形状を表す形状パラメータを推定する。推定された形状パラメータが表す特定の物体の形状を事前知識として、被写体の領域を、入力画像から抽出する。これにより、計算量の増大を抑制して、被写体の領域を精度よく抽出することができる。

【0148】
[第4の実施の形態]
次に、第4の実施の形態について説明する。なお、第4の実施の形態に係る統計的形状モデル生成装置及び画像処理装置の構成は、第1の実施の形態と同様の構成となるため、同一符号を付して説明を省略する。

【0149】
第4の実施の形態では、目的関数を拡張させる点が、第1~第3の実施の形態と異なる。

【0150】
上記第1の実施の形態では、上記式(4)を目的関数として用い、上記式(4)におけるF(I,y)及びB(I,y)は、上記式(5)、(6)で定義される場合について説明した。

【0151】
第4の実施の形態では、上記式(4)に示す目的関数は、特定の物体の形状の尤もらしさとして、形状パラメータが表す特定の物体の形状に対する画素における形状ラベルの値に関して単調に変化する単調関数を含むように定義し、目的関数を拡張する。

【0152】
具体的には、第4の実施の形態では、上記式(4)におけるF(I,y)とB(I,y)とを、以下の式(18)、(19)に示すように定義し、目的関数を拡張する。

【0153】
【数15】
JP2016027840A1_000018t.gif

【0154】
図16に、目的関数の下界が単調性を満たすための十分条件を表す図を示す。この十分条件は,h(y),h(y)(∀p∈P)がそれぞれy(∀q∈P)に関して単調減少,単調増加する単調関数であることである。図16に示すように、ラベル対{y,y’}∈L|P|

【0155】
【数16】
JP2016027840A1_000019t.gif

【0156】
であるとき、以下の式(20)に示す関係が成り立つ。

【0157】
【数17】
JP2016027840A1_000020t.gif

【0158】
また、図17に、上記式(20)の関係を満たすh(y)とh(y)の一例を示す。図17に示す例は、以下の式(21)に示す距離関数によってh(y)とh(y)とを定義したものである。

【0159】
【数18】
JP2016027840A1_000021t.gif

【0160】
上記式(21)に示す符号付き距離関数によって単調関数h(y)とh(y)とを定義し、形状の輪郭付近でのペナルティを小さくすることで、事前形状の誤りによるセグメンテーションの誤差を低減することができる。

【0161】
以上説明したように、第4の実施の形態に係る画像処理装置によれば、目的関数における特定の物体の形状の尤もらしさとして、形状パラメータが表す特定の物体の形状に対する画素における形状ラベルの値に関して単調に変化する単調関数を含むように定義する。これにより、被写体の領域を精度よく抽出することができる。

【0162】
また、第4の実施の形態では、特定の物体の形状の正解の輪郭でエネルギーが最小になるような単調関数h(y)とh(y)とを目的関数に導入することで、より優れた目的関数が利用可能となる。

【0163】
[第5の実施の形態]
次に、第5の実施の形態について説明する。なお、第5の実施の形態に係る統計的形状モデル生成装置及び画像処理装置の構成は、第1の実施の形態と同様の構成となるため、同一符号を付して説明を省略する。

【0164】
第5の実施の形態では、固有空間上のパラメータαを形状に写像する関数gを拡張させる点が、第1~第4の実施の形態と異なる。

【0165】
上記第1の実施の形態において、上記式(2)に写像の関数gは、上記式(3)のHeaviside function H(・)によって定義される場合について説明した。

【0166】
第5の実施の形態では、形状パラメータを特定の物体の形状に写像する関数gは、予め定められた関数fを含むように定義し、写像する関数gを拡張する。

【0167】
【数19】
JP2016027840A1_000022t.gif

【0168】
また、写像する関数gを以下の式(23)のように拡張した結果、レベルセット関数以外の形状表現も利用可能となる(例えば、Kilian M. Pohl et al.,“Logarithm Odds Maps for Shape Representation”, Proc. of Medical Image Computing and Computer-Assisted Intervention, Vol.4191, pp.955-963, 2006 に記載されているLogOddsなど)。

【0169】
【数20】
JP2016027840A1_000023t.gif

【0170】
ここでfは予め定められた単調関数を表し、tは予め定められた閾値を表す。

【0171】
例えば、統計的形状モデルが表す可能な形状表現の一例として、LogOddsを用いる場合について説明する。形状を表すクラスが2である場合、統計的形状モデルとしてLogOddsを用いるときは、以下の式で表される。

【0172】
【数21】
JP2016027840A1_000024t.gif

【0173】
また、図18に、統計的形状モデルとしてLogOddsを用いたとき例を示す。図18に示すように、LogOddsを用いることにより、形状表現が拡張される。

【0174】
以上説明したように、第5の実施の形態に係る統計的形状モデル生成装置によれば、写像する関数gの拡張(予め定められた関数fと閾値tとの導入)により、さらに高度な統計的形状モデルが利用可能となる。

【0175】
なお、本発明は、上述した例に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

【0176】
本実施の形態は、基礎的であることから、その応用範囲は広いと考えられる。具体的には、形状などが統計的に変動する図形を対象とした画像認識処理に幅広く利用可能である。例えば、他の臓器や他の医用画像(例えばMRやPETなど)はもちろん、それ以外に、顔画像からの顔の認識や、一般の風景画像から特定の図形を認識する処理などに利用可能である。本実施の形態では、膵臓を対象としたが、他の臓器として例えば脾臓を対象とすることもできる。

【0177】
また、第1の実施の形態では、統計的形状モデルとしてLSDMを用いる場合を例に説明したが、これに限定されるものではない。LSDMでなく、他の統計的形状モデルの場合にも利用可能である。統計的形状モデルは、上記式(1)のような線形の関数(厳密には画素についての関数φが線形の関数)で表現できていれば、他の統計的形状モデルに本実施の形態を適用することができる。目的関数や最適化手法によっては、実施の形態で示した上記式(1)のような線形の関数で表現される統計的形状モデル以外の統計的形状モデルに対しても、本実施の形態を適用することができる。
また、上記第5の実施の形態に示すように、写像の関数gを拡張させることで、様々な統計的形状モデルを用いることができる。

【0178】
また、第1の実施の形態では、固有空間の次元は2(d=2)である場合を例に説明したが、これに限定されるものではない。例えば、より高次元の固有空間を対象としてもよい。例えば、第2及び第3の実施の形態で示したように、近似解法を用いることによって、固有空間の次元を増加させることができる。

【0179】
上記実施の形態において、固有空間の次元を3以上に増加させた場合、凸多角形は凸多胞体に対応し、凸多角形の面積は、凸多胞体の体積に対応する。
従って固有空間の次元を3以上に増加させた場合、セグメンテーション部114は、固有空間上の凸多胞体を繰り返し分割して最適な形状パラメータを示す点を含む凸多胞体を探索する探索アルゴリズムに従って、上記式(4)に示す目的関数を最適化するように、入力画像が表す被写体の形状を表す形状パラメータを推定する。
また、セグメンテーション部114は、探索アルゴリズムにおいて、凸多胞体に含まれる形状集合に対する目的関数の下界を計算することにより、最適な形状パラメータを示す点を含む固有空間上の凸多胞体を探索して、形状パラメータを推定する。
また、セグメンテーション部114は、分割後の2つの凸多胞体の体積が対応するように、固有空間上の凸多胞体を分割し、かつ、最適な形状パラメータを示す点を含む固有空間上の凸多胞体を探索することを繰り返して、形状パラメータを推定する。

【0180】
また、実施の形態では、上記式(17)に従って選択された集合Qに含まれる画素kによって固有空間を分割する場合を例に説明したが、これに限定されるものではなく、他の任意の分割法で作成した凸多角形に対しても適用可能である。例えば、第2の実施の形態で示したように、格子状など、規則的に固有空間を分割した凸多角形にも適用可能である。
この場合、分割の細かさによってはTightness(最適性保証の3つの条件の内の一つ。)が成立しないことがある。

【0181】
ここで、Tightnessについて簡単に説明する。

【0182】
上記式(14)を、以下の式(24)のように示す。

【0183】
【数22】
JP2016027840A1_000025t.gif

【0184】
この場合、すべての葉ノードHについて以下の式(25)が成立する場合には、Tightnessが成立している。

【0185】
【数23】
JP2016027840A1_000026t.gif

【0186】
分割の細かさが粗い場合には、上記式(25)が成立せずに、Tightnessが成立しないこともあり、近似解しか得られないことがあるが、固有空間の次元dを大きくした場合の計算コストのさらなる削減が可能であるというメリットがある。注意すべき点は、最適性をやや犠牲にすることで、計算コストと計算精度のバランスをユーザが選ぶことが可能になる点であり、本実施の形態の提案法の拡張に関する重要な点である。なお、このとき、最適性が保証されない点では上記米国特許第8249349号明細書などと同様であるが、上記米国特許第8249349号明細書と比べて、前処理が一切不要である点に注意されたい。

【0187】
また、最小化の目的関数やアルゴリズムについても、今回は上記実施の形態で示した目的関数、および、branch and bound法とグラフカット法とを用いたが、最小化が可能であれば、今回示したものに限らず適用可能である。

【0188】
また、上記第2又は第3の実施の形態に係る近似解法を実行した後に、上記第1の実施の形態に係る厳密解法を実行してもよい。

【0189】
また、図2に示したコンピュータ構成において、実施の形態に係る各処理部の機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより、各構成による処理が実行されてもよいし、図示されていない通信機能を用いて、当該プログラムを読み込ませることでもよい。

【0190】
なお、コンピュータ読み取り可能な記録媒体とは、フレキシブルディスク、光磁気ディスク、ROM、CD-ROM等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。

【0191】
また、上記プログラムは、このプログラムを記憶装置等に格納したコンピュータシステムから、伝送媒体を介して、あるいは、伝送媒体中の伝送波により他のコンピュータシステムに伝送されてもよい。ここで、プログラムを伝送する「伝送媒体」は、インターネット等のネットワーク(通信網)や電話回線等の通信回線(通信線)のように情報を伝送する機能を有する媒体のことをいう。

【0192】
また、上記プログラムは、前述した機能の一部を実現するためのものであってもよい。さらに、前述した機能を、コンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるもの、いわゆる差分ファイル(差分プログラム)であってもよい。

【0193】
また、本実施の形態では、上記図1に示す統計的形状モデル生成装置10及び上記図7に示す画像処理装置100の各処理部を、プログラムにより各機能の実行が可能なコンピュータで構成するものとしているが、論理素子回路からなるハードウェア構成とすることでも良い。

【0194】
また、上記図2に示す統計的形状モデル生成装置10及び画像処理装置100のコンピュータ構成に関しても適宜にその構成内容を変更しても良い。

【0195】
このように、実施する形態例を、図面を参照して詳述してきたが、具体的な構成はこの実施の形態例に限られるものではなく、要旨を逸脱しない範囲の設計等も含まれる。

【0196】
また、実施の形態のプログラムは、記憶媒体に格納して提供するようにしてもよい。

【0197】
実施の形態のコンピュータ可読媒体は、特定の物体である被写体を表す入力画像から、前記被写体の領域を抽出するプログラムであって、コンピュータを、前記入力画像を受け付ける受付手段、及び(A)前記被写体を表す学習用の複数の画像であって、前記被写体の領域が予め求められた前記複数の画像に基づいて予め計算された固有ベクトルを基底とする固有空間であって、かつ、(B)固有空間上の点が、前記特定の物体の形状の統計的変動を表す統計的形状モデルの形状パラメータを示す固有空間において、前記入力画像に基づいて、前記固有空間上の点が示す前記形状パラメータが表す前記特定の物体の形状の尤もらしさと前記入力画像中の隣接画素間の画素値の差とに応じた値を表す予め定められた目的関数を最適化するように、前記入力画像が表す前記被写体の形状を表す前記形状パラメータを推定し、前記推定された前記形状パラメータが表す前記特定の物体の形状を事前知識として、前記被写体の領域を、前記入力画像から抽出するセグメンテーション手段として機能させるためのプログラムを記憶したコンピュータ可読媒体である。

【0198】
日本出願2014-169911の開示はその全体が参照により本明細書に取り込まれる。

【0199】
本明細書に記載された全ての文献、特許出願、及び技術規格は、個々の文献、特許出願、及び技術規格が参照により取り込まれることが具体的かつ個々に記載された場合と同程度に、本明細書中に参照により取り込まれる。
図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13
【図15】
14
【図16】
15
【図17】
16
【図18】
17
【図19】
18