JP3649328B2 - 画像領域抽出方法および装置 - Google Patents
画像領域抽出方法および装置 Download PDFInfo
- Publication number
- JP3649328B2 JP3649328B2 JP2002108608A JP2002108608A JP3649328B2 JP 3649328 B2 JP3649328 B2 JP 3649328B2 JP 2002108608 A JP2002108608 A JP 2002108608A JP 2002108608 A JP2002108608 A JP 2002108608A JP 3649328 B2 JP3649328 B2 JP 3649328B2
- Authority
- JP
- Japan
- Prior art keywords
- class
- coarse
- grained
- image
- subspace
- 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.)
- Expired - Fee Related
Links
Images
Classifications
-
- 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/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
- G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
- G06F18/2415—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on parametric or probabilistic models, e.g. based on likelihood ratio or false acceptance rate versus a false rejection rate
-
- 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/143—Segmentation; Edge detection involving probabilistic approaches, e.g. Markov random field [MRF] modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/60—Type of objects
- G06V20/69—Microscopic objects, e.g. biological cells or cellular parts
- G06V20/698—Matching; Classification
-
- 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/10024—Color image
-
- 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/10056—Microscopic image
-
- 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/30024—Cell structures in vitro; Tissue sections in vitro
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Probability & Statistics with Applications (AREA)
- Data Mining & Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Evolutionary Computation (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Artificial Intelligence (AREA)
- Health & Medical Sciences (AREA)
- Biomedical Technology (AREA)
- General Health & Medical Sciences (AREA)
- Molecular Biology (AREA)
- Multimedia (AREA)
- Software Systems (AREA)
- Image Analysis (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、領域抽出の対象像を含む画像データに対し、画像のピクセルが持つ属性に基づいて、対象となる物体の領域を抽出する技術に関し、特に、ＭＲＩ装置やＣＴスキャン装置などで撮影された画像データ、あるいは生物顕微鏡で観測された生体組織の断面図から、画像を特定の臓器や腫瘍、細胞の核、腺腔などの領域を抽出する画像領域抽出方法および装置に関する。
【０００２】
【従来の技術】
従来、画像から特定の物体の像を抽出するには、次のような方法が知られている（例えば、特開２００１−９２９８０号公報；以下文献１とする。）
文献１にでは、輪郭抽出の対象となる物体の像を含む画像データに対し、前記画像の各点が持つ属性に基づいて、画像各点が各領域に属する確率である領域帰属確率を用いて、前記各点が属する領域を分離し、前記領域間の境界を輪郭として抽出する構成としている。このような構成とすることにより、領域分類の閾値を陽（explicit）に設定することなく、自動的に各領域の輪郭を抽出することが可能となり、輪郭領域抽出の高速化を可能にしている。
【０００３】
【発明が解決しようとする課題】
上記文献１における輪郭領域抽出は、画像データから、画像上各点が領域の各々に属する領域帰属確率の期待値を算出し、領域帰属確率と領域パラメータから求められる混合確率分布に基づいて評価関数を算出し、領域帰属確率に基づいて各領域を分離し、分離された領域に基づいて輪郭を抽出しているため、評価関数を求める際、すべての画像点（ピクセル）に関して和を取る必要があり、最適なパラメータを求めるためにはこの評価関数の算出を何度も繰り返すことになる。
【０００４】
そのため、従来の方法においては，画像のサイズが大きくなればなるほど、領域抽出に要する時間が膨大になるという問題がある。
【０００５】
本発明は上述した問題点を鑑み、画像から、閾値を陽に設定することなく、更なる高速かつ高精度に対象領域を抽出する画像領域抽出方法および装置を提供することにある。
【０００６】
【課題を解決するための手段】
本発明第１記載の画像領域抽出方法は、
領域抽出の対象となる画像データに対し、画像上の各ピクセルを複数のクラスに分割し、所望の領域を抽出する画像領域抽出方法において、
前記画像の前記各ピクセルが取りうる属性値全体から構成されるデータ空間を、所与の解像度で部分空間に分解し、該各部分空間に属性値を取るピクセルの番号と、該ピクセルの属性値の平均値と、該ピクセルの個数とを保持し、粗視化データ空間を構成する第１の工程と、
前記各部分空間のピクセルの個数を前記画像に含まれる全ピクセルの個数で割り、前記粗視化データ空間上の粗視化経験確率分布を算出する第２の工程と、
前記各クラスの属性を規定するクラスパラメータ、前記クラスの数、および前記クラスの混合比率を初期化する第３の工程と、
前記クラスパラメータに基づいて、前記各部分空間の前記平均値における条件付き確率を算出し、該条件付き確率のデータ空間における分布である粗視化条件付き確率分布を算出する第４の工程と、
前記画像を構成する各ピクセルが各前記クラスに帰属する確率であるクラス帰属確率を、前記混合比率に前記粗視化条件付き確率分布を掛けることにより算出する第５の工程と、
粗視化対数尤度を評価関数として用い、該評価関数を増加させるように前記クラスパラメータ、および前記クラスの混合比率を更新する第６の工程と、
前記粗視化条件付き確率分布と、前記粗視化経験確率分布と、前記第６の工程で更新した前記クラスパラメータおよび前記クラスの混合比率とを用いて、前記評価関数を算出する第７の工程と、
前記評価関数が所与の終了条件を満たしているかどうかを調べる第８の工程と、
前記評価関数が前記所与の終了条件を満たした後、前記クラスパラメータおよび前記混合比率を保持し、前記クラス帰属確率に基づいて、前記ピクセルが属するクラスを決定し、前記所望の領域を抽出する第９の工程とを有し、
前記評価関数が前記所与の条件を満たすまで、前記第４、前記第５、前記第６、前記第７、および前記第８の工程を繰り返すことを特徴とする。
【０００８】
本発明第２に記載の画像領域抽出方法は、
前記第１に記載の画像領域抽出方法において、粗視化の解像度が元の解像度であるかを調べる第１０の工程と、前記部分空間の解像度を元の解像度へ復元する第１１の工程とを有し、
前記第３の工程において、前記第９の工程で保持されているクラスパラメータおよび前記クラスの混合比率を初期値として用い、所与の条件が満たされるまで、前記第４、前記第５、前記第６、前記第７、および前記第８の工程を繰り返すことを特徴とする。
【０００９】
本発明第３記載の画像領域抽出方法は、
前記第１及至第２のいずれかに記載の画像領域抽出方法における前記第９の工程において、推定された前記混合比率に前記画像を構成する全ピクセルの個数を掛けて前記各クラスに属するピクセルの個数を算出し、前記クラス帰属確率が高い順から前記ピクセルを選択し、前記各クラスに属するピクセルを決定することを特徴とする。
【００１０】
本発明第４記載の画像領域抽出方法は、
前記第１及至第３のいずれかに記載の画像領域抽出方法における前記第６の工程において、前記評価関数としてＡＩＣ（Akaike Information Criterion）を用い、前記第７の工程において評価関数を減少させるようパラメータを変更することを特徴とする。
【００１１】
本発明第５記載の画像領域抽出方法は、
前記第１及至第４のいずれかに記載の画像領域抽出方法における前記第６の工程において、前記評価関数としてＭＤＬ（Minimum Description Length）を用い、前記第７の工程において評価関数を減少させるようパラメータを変更することを特徴とする。
【００１２】
本発明第６記載の画像領域抽出方法は、
前記第１及至第４のいずれかに記載の画像領域抽出方法における前記第６の工程において、前記評価関数としてストラクチュラル・リスク（Structural Risk）を用い、前記第７の工程において該評価関数を減少させるようパラメータを変更することを特徴とする。
【００１３】
本発明第７記載の画像領域抽出方法は、
前記第１及至第６のいずれかに記載の画像領域抽出方法における前記第３の工程において、
前記各部分空間同士が近傍にあるかどうかを定義する近傍半径と、前記クラスの数を設定し（ステップＢ１）、
前記各部分空間ごとに、各部分空間の代表値を設定し（ステップＢ２）、
分類対象部分空間の集合を設定し（ステップＢ３）、
前記分類対象部分空間の集合の中から、前記粗視化経験確率の最も高い部分空間を選択し（ステップＢ４）、
前記粗視化経験確率の最も高い部分空間の代表値との距離が、前記近傍半径以内である代表値を持つ部分空間全てを、近傍集合として選択し（ステップＢ５）、
すでに分類が完了したクラスに含まれる部分空間の代表値と、前記近傍集合に含まれる部分空間の代表値との間の最短距離が、前記近傍半径よりも大きいかどうかを調べ（ステップＢ６）、
もし、すでに分類が完了したクラスに含まれる部分空間の代表値と、前記近傍集合に含まれる部分空間の代表値との間の最短距離が、前記近傍半径よりも大きければ、前記近傍集合を新規クラスとし、前記分類対象部分空間から前記近傍集合を削除し、ステップＢ４以降を繰り返し（ステップＢ７）、
そうでなければ、前記近傍集合をすでに分類済みのクラスに追加し、前記分類対象部分空間から、前記近傍集合を削除し（ステップＢ８）、
前記分類対象部分空間が空集合であるかどうかを調べ（ステップＢ９）、
もし、前記分類対象部分空間が空集合でなければ、ステップＢ４以下を繰り返し、
前記分類対象部分空間が空集合であれば、分類が完了したクラスの数が所与の数以上あるかどうかを調べ（ステップＢ１０）、
もし分類が完了したクラスの数が所与の数よりも少なければ、前記近傍半径を縮小し（ステップＢ１１）、ステップＢ３以下を繰り返し、
もし、前記分類対象部分空間が空集合であり、かつ分類済みのクラスの数が所与の数よりも多ければ、各クラス内で前記クラスパラメータを算出し、これをクラスパラメータの初期値とし、また各クラス内に含まれる部分空間の数の比率を前記混合比率の初期値とする（ステップＢ１２）、ことを特徴とする。
【００１４】
本発明第８記載の画像領域抽出装置は、
領域抽出の対象となる像を含む画像データに対し、画像上の各ピクセルを複数のクラスに分類し、所望の領域を抽出する画像領域抽出装置において、
画像データを読み込む入力装置と、
前記画像の前記各ピクセルが取りうる属性値全体から構成されるデータ空間を、所与の解像度で部分空間に分解し、該各部分空間に属性値を取るピクセルの番号と、該ピクセルの属性値の平均と、該ピクセルの個数を保持し、粗視化データ空間を構成する領域粗視化装置と、
前記各部分空間のピクセルの個数を前記画像に含まれる全ピクセルの個数で割り、前記粗視化データ空間上の粗視化経験分布を算出する粗視化経験確率分布算出装置と、
前記各クラスの属性を規定するクラスパラメータ、前記クラスの数、および前記クラスの混合比率を初期化し、前記各クラスの属性を規定するクラスパラメータから、前記クラスが指定されたもとでの条件付き確率分布を算出し、前記クラスが指定されたもとでの条件付き確率分布を、前記各部分空間内で平均化することによって粗視化条件付き確率分布を算出する粗視化条件付き確率分布装置と、
前記画像を構成する各ピクセルが前記各クラスに帰属する確率であるクラス帰属確率を、前記粗視化条件付き確率分布から算出するクラス帰属確率算出装置と、
評価関数として粗視化対数尤度を用いるときには、該評価関数を増加させるように、また評価関数としてＡＩＣ、ＭＤＬ、ストラクチュラル・リスクのいずれかを用いる場合は該評価関数を減少させるように前記パラメータを更新するパラメータ更新装置と、
前記粗視化対数尤度、前記ＡＩＣ、前記ＭＤＬ、前記ストラクチュラル・リスクのいずれかを用いて前記評価関数を算出する評価関数算出装置と、
前記評価関数が所与の終了条件を満たしているかどうかを調べ、前記評価関数が前記所与の終了条件を満たした後、前記パラメータを保持し、前記クラス帰属確率に基づいて、前記各ピクセルが属するクラスを決定し、領域を抽出する領域抽出装置と、
抽出した領域を出力する出力装置と、
を有することを特徴とする。
【００１５】
本発明第９記載の画像領域抽出装置は、前記第８記載の画像領域抽出装置において、前記評価関数が前記所与の終了条件を満たしていることが確認されたのち、粗視化の解像度が元の解像度であるかを調べ、前記データ空間の解像度を元の解像度へ復元する解像度復元装置を有することを特徴とする。
【００１６】
また、本発明は、本発明における画像領域抽出処理をコンピュータに実行させるための画像領域抽出プログラムをその特徴とする。
【００１７】
本発明では文献１記載の発明と同様に、推定したクラス帰属確率に基づいて領域を抽出している。この特徴により、閾値を自動的に設定できるという利点がある。しかしながら、このクラス帰属確率を推定するためには、パラメータの更新を何度も（通常５０回程度）繰り返し行う必要があり、計算に時間を要するという問題がある。
【００１８】
そこで、本発明では画像データの発生源のモデルとして粗視化確率分布を用い、それに伴ってデータ空間も粗視化する。これにより、通常は条件付き確率分布やクラス帰属確率を数十万個のデータについて計算する必要があるところを、数千個の部分空間についての計算に置き換えることができる。この特徴により、領域抽出に要する時間を大幅に短縮することが可能となる。
【００１９】
推定されるクラスパラメータは、粗視化によって精度が劣ることになるが、本発明ではクラスパラメータを直接使うのではなく、パラメータから計算したクラス帰属確率に基づいて領域抽出を行うため、後に示すように、粗視化による誤差の影響をうけることなく領域抽出を行うことが可能となっている。
【００２０】
以下、上記本発明の作用の有効性について更に詳細に説明する。
【００２１】
本発明の第１に記載の画像領域抽出方法においては、画像を構成する各ピクセルの持つ属性値を確率変数と見なし、推定したピクセル値の確率分布に基づいて所望の領域を抽出する。属性値としては、例えばモノクロ画像ならば輝度の値、カラー画像ならば赤（Ｒ）、緑（Ｇ）、青（Ｂ）の色要素の強度などを用いることができる。
【００２２】
所望の領域を抽出するためには、各ピクセルの属性値に基づいて、類似の属性をもった複数のグループに各ピクセルを分類する必要がある。本明細書においては、類似の属性を持ったピクセルの集合をクラスと呼ぶ。各クラスは、例えばそのクラスに属するピクセルの属性値の平均値や分散などによって特徴づけられるものとする。以下、これらクラスの特徴をそのクラスの「クラスパラメータ」と呼び、φi（i=1,,,k）と書くこととする。ｋはクラスの数である。このような定式化のもとで、ｊ番目のピクセルが、ｘjという属性値を取る確率は、次の混合分布で表現できる。
【００２３】
【数１】
ここで、f(xj|φi)はデータがｉ番目のクラスから発生していると仮定したときの条件付確率分布、ｗiは、各クラスの混合比で、
【００２４】
【数２】
を満たす。
【００２５】
ｘjは、例えばモノクロ画像ならば、0から255などの整数値をとり、カラー画像ならばＲＧＢの色要素の値を成分とした３次元ベクトル（xj1,xj2,xj3）で表される。ここで各ｘjl（l=1,2,3）は、例えば0から255の整数値を取る。
【００２６】
混合比ｗiは、異なるクラスに属する領域の面積比を表す。例えば、平均輝度200、輝度の標準偏差が20という特徴で表される明るい画像領域（クラス１とする）と、平均輝度50、輝度の標準偏差が10という特徴で表される暗い領域（クラス２とする）の、二つの領域から構成されるモノクロ画像があるとする。また、その画像の面積の７割を明るい領域が、３割を暗い領域が占めているとする。この場合、クラスの数はｋ=2、クラスパラメータはφ1=(200,20),φ2=(50,10)であり、この画像の混合分布は、
p(xj)=0.7(xj|200,20)+0.3(xj|50,10) ……(２）
などと表現できる。
【００２７】
以下、クラス混合比ｗiとクラスパラメータφiをまとめてθiと書くことにする。以下において単に「パラメータ」というときは、このθiを意味することとする。
【００２８】
本発明第１記載の画像領域抽出方法では、以下のように定義される平均対数尤度
【００２９】
【数３】
を最大にするパラメータを推定し、推定したパラメータの情報を利用して領域抽出を行う。ここでｎは画像に含まれるピクセル数である。このような統計手法は最尤法と呼ばれている。
【００３０】
しかしながら、平均対数尤度Ｌを最大にするパラメータを推定することは、一般に困難である。そこで、この平均対数尤度Ｌの代わりに次の量
【００３１】
【数４】
で表される完全対数尤度の期待値Ｑを用いてパラメータを推定することができる。ここで
【００３２】
【数５】
は、ｊ番目のピクセルがｉ番目のクラスに帰属する確率である。本発明ではこれをクラス帰属確率と呼ぶ。Ｑが増大するようにパラメータを更新すれば、上述の平均対数尤度Ｌも必ず増大することは、数学的に証明されており、この証明は、例えば A. P. Dempster, N.M.Laird, and D.B.Rubin, Maximum Likelihood From Incomplete Data via The EM Algorithm, J. Roy. Stat. Soc. Vol.30, pp.205 248, 1977（以下文献２とする）に詳しい。
【００３３】
本発明において、推定されたパラメータから実際に領域を抽出するには、以下のように行う。まず適当な初期パラメータから開始し、（５）式で与えられるクラス帰属確率を求める。次にＱを増加させるようにパラメータｗ，φを更新し、改めてＱを算出する。最終的にＱが増加しなくなるまでこの手順を繰り返す。
【００３４】
ｋ個のクラスのうち、例えばｉ番目のクラスに属する点を抽出するには、ｉ番目のクラスへの帰属確率の値を調べ、その確率値が一定値以上のピクセルをそのクラスに属するものとして分類する。すべてのピクセルを分類し終わった時点で、ｋ個のクラスの中から所望の属性を持ったクラスを選択し、そのクラスに属するピクセルを画像から抽出すれば、所望の領域を自動的に抽出することができる。
【００３５】
本発明では、このＱの最大化を高速に行うために、粗視化確率分布を導入する。粗視化確率分布は、データが取りうる値全体からなる空間（以下、データ空間と呼ぶ）を互いに交わらないＮ個の部分空間に分解し、各部分空間に確率を割り当てることによって構成される。
【００３６】
具体的には、ｊ番目の部分空間における粗視化された条件付き確率分布を
【００３７】
【数６】
で定義し、粗視化条件付き確率分布を
【００３８】
【数７】
とする。ここで、Ａjはｊ番目の部分空間である。Ｄをデータ空間全体とすると
【００３９】
【数８】
を満たす。またＩＡ(x)は、データ値が部分空間Ａに含まれるときに１、それ以外の時に０となる指示関数で、m(A)=∫Adx はＡの測度である（データ空間が２次元のときにはＡの面積を、３次元空間のときは体積を表す。）
このように定義された粗視化条件付き確率分布を用いると、Ｑは、
【００４０】
【数９】
と書き換えることができる。ここで
【００４１】
【数１０】
は粗視化された経験確率分布、
【００４２】
【数１１】
は粗視化されたクラス帰属確率である。式（９）で与えられる粗視化された完全対数尤度の期待値を最大化することによって、次の粗視化平均対数尤度
【００４３】
【数１２】
を最大化することが可能となる。
【００４４】
元のＱと比較すると、式（４）で与えられるＱが全データに関して和をとっているのに対し、式（９）で与えられる粗視化された完全対数尤度は、部分空間に関してのみ和を取っている。後に示すように、この改善によって本発明では計算量が大幅に軽減できる。例えば、512×512ピクセルの画像の場合、２６万個以上のデータについて和を取る必要があるが、本発明の粗視化分布を用いた場合、部分空間に関する和を１０００個程度にまで軽減でき、高速な推定が可能となる。
【００４５】
上記の方法において各部分空間における粗視化された確率値は、その部分空間に含まれるピクセルの属性値の平均値における確率値で近似する。
【００４６】
【数１３】
ここで
【００４７】
【数１４】
は、ｊ部分空間Ａjに含まれるピクセルの属性値の平均値である。この近似によって、部分空間内における積分（あるいは総和）操作を省略することができ、さらなる計算量の軽減が可能となる。
【００４８】
本発明第２に記載の画像領域抽出方法では、粗視化確率分布を用いて推定したパラメータを初期値とし、元の解像度で再度推定を行う。この場合、粗視化確率分布を用いてすでにほぼ最適なパラメータが得られているため、最初から元の解像度で推定を行う場合に比較して、パラメータの逐次更新の回数がはるかに少なくてすむため、高精度な推定を高速に行うことが可能となる。
【００４９】
本発明第３に記載の画像領域抽出方法では、領域を抽出する際、推定された混合比ｗiに全ピクセル数ｎを掛け、ｉ番目の領域に属するピクセル数を推定する。そして領域帰属確率の高い順に上位ｎi個をこの領域に属するピクセルとして抽出する。この方法により、どの値までの確率値をその領域に属しているものと見なすか、という閾値を自動的に決定することができる。
【００５０】
本発明第４、第５および第６に記載の画像領域抽出方法においては、評価関数として、それぞれＡＩＣ、ＭＤＬ、ストラクチュラル・リスクを用い、その値がもっとも小さい結果を与えるモデルを選択する。最尤推定においては一般にパラメータの数が多いほど、評価値としてはよい値を与えるため、本来のパラメータ数が決定できないという問題がある。
【００５１】
これに対し、本発明第４、第５および第６に記載の評価関数は、過剰なパラメータ数を用いたとき、評価関数の値が逆に増加するという特性を持っている。これにより、最適なパラメータ数を推定することが可能となる。領域抽出へ適用する場合は、画像が何種類の領域から形成されているか、その適切な領域数を推定することができる。
【００５２】
【発明の実施の形態】
次に、本発明の好ましい実施の形態について、図面を参照しながら詳細に説明する。なお、以下の説明において用いる記号は、特に言及しない限りこれまで用いてきた記号の説明に従うものである。また、以下の実施例では、染色された細胞のカラー顕微鏡写真から、細胞の各組織（核、腺腔）を抽出する場合を例に取って説明を行うが、画像として記録されている像であるならば、以下に説明する方法と同様な方法で実行できる。
【００５３】
図１は本発明の第１の実施形態の画像領域抽出方法における処理手順を示すフローチャートである。また、図２は本実施形態の画像領域抽出装置の構成を示すブロック図である。
【００５４】
本発明の画像領域抽出装置は、画像データを読み込む入力装置１と、データ空間粗視化装置２と、粗視化経験確率分布算出装置３と、粗視化条件付き確率分布算出装置４と、クラス帰属確率算出装置５と、パラメータ更新装置６と、評価関数算出装置７と、領域抽出装置８と、出力装置９から構成されている。
【００５５】
以下、本発明の画像領域抽出装置の処理について、図１および図２を用いて説明する。
【００５６】
入力装置１は、生物顕微鏡などで観測された細胞のカラー画像を入力する（ステップＡ１）。この入力装置は、例えば、画像スキャナやデジタルカメラなどを用いて構成することができる。あるいは、生物顕微鏡をコンピュータに接続し、ネットワークを通して画像を直接入力することも可能である。入力装置１は、読み込んだデータをデータ空間粗視化装置２に送る。
【００５７】
データ空間粗視化装置２は、本発明第１記載の第１の工程に基づいて、粗視化データ空間を構成する（ステップＡ２）。ここで、データ空間とは、画像の各ピクセルが取りうる属性値全体の集合を意味する。例えば標準的グレースケールの画像において、各ピクセルの属性は輝度によって表すことができ、輝度は通常0から255の整数値で表現される１次元空間である。また、標準的なカラー画像の場合ならば、Ｒ、Ｇ、Ｂ各色要素に対し、通常0から255の整数値で表現される３次元空間となる。
【００５８】
部分空間を構成するには、例えば解像度を８とすると、ＲＧＢの各値を８ずつに区切り、８×８×８の立方体を１つの部分空間として構成することができる。解像度は３次元の軸で同一で無くともよい。例えばＲ軸をｈ1ずつ、Ｇ軸をｈ2ずつ、Ｂ軸をｈ3ずつに区切って粗視化してもよく、各部分空間が重なりを持たず、かつ全データ空間をもれなく覆うことができればよい。以下において、解像度ｈで粗視化すると表現する場合は、ｈ×ｈ×ｈでＲＧＢの各値を区切ることを意味する。
【００５９】
データ空間粗視化装置２は、また、各部分空間に値を取るピクセルの集合と、それらピクセルの属性値の平均と、ピクセルの個数を保持する。
【００６０】
図６に解像度ｈ1×ｈ2×ｈ3で粗視化したときの粗視化データ空間１１、および部分空間１２を示す。例えばｊ番目のピクセルの属性値ｘjが、Ｒ=128、Ｇ=255、Ｂ=18なる値の時、このピクセルは、（16,31,2）なる指標で指定される部分空間に値を取る。このことを以下では、ｘjがこの部分空間に含まれる、と表現する。
【００６１】
データ空間粗視化装置２は、粗視化データ空間を構成した後、各部分空間に含まれるピクセルの番号と、ピクセル数を粗視化経験確率分布算出装置３に送る。粗視化経験確率分布算出装置３では、本発明第１記載の第２の工程に基づいて、粗視化経験確率分布を算出する（ステップＡ３）。ここで粗視化経験確率分布とは、粗視化データ空間の各部分空間に含まれるピクセルの個数を、全ピクセル数で割った値から算出された確率分布を表す。
【００６２】
粗視化経験確率分布算出装置３は、粗視化経験確率分布を粗視化条件付き確率分布算出装置４に送る。粗視化条件付き確率分布算出装置４では、本発明第１記載の第３の工程に基づいて、クラスパラメータを初期化する（ステップＡ４）。また、粗視化条件付き確率分布算出装置４では、本発明第１記載の第４の工程に基づいて、粗視化条件付き確率分布を算出する（ステップＡ５）。
【００６３】
粗視化条件付き確率分布を具体的に算出するには、以下のようにする。ここでは、ｊ番目のピクセル値がｉ番目のクラスから発生しているとした条件のもとでの条件付き確率が、以下のような多次元正規分布
【００６４】
【数１５】
で与えられているとする。ここで、ｘは、ＲＧＢ各色の値からなる３次元ベクトル、μiはｉ番目のクラスの平均色を表す３次元ベクトル、Σiはｉ番目のクラスの共分散行列であり、|Σi|,Σi -1はそれぞれ行列Σiの行列式，逆行列を表す。また、(x-μi)Ｔは転置を表す。
【００６５】
粗視化条件付き確率分布算出装置４では、この条件付き確率を式（６）で与えられる式で計算する。このとき、各部分空間の測度ｍ(Ａj)は、各部分空間の体積となる。例えば、解像度８で一様に粗視化した場合、8×8×8=512となる。
【００６６】
粗視化条件付き確率分布算出装置４では、本発明第１記載の方法のように、粗視化条件付き確率分布を式（１３）によって近似することも可能である。この方法により、パラメータを更新する度に式（６）で与えられる演算をする必要がなくなり、計算量を大幅に削減できる効果がある。
【００６７】
また、粗視化条件付き確率分布算出装置４では、本発明第７記載の方法のように、粗視化経験分布に基づいてパラメータの初期値を決定することも可能である。本発明第７記載の方法では、各部分空間を大雑把にクラス分類し、分類の結果得られた各クラス内で平均値や分散値を求め、これらの値をパラメータ推定の初期値とする。
【００６８】
図５は、粗視化経験確率分布に基づいて、パラメータの初期値を決定する処理手順を示したフローチャートである。以下、図５を参照しながらパラメータの初期値設定法について説明する。
【００６９】
まず、ステップＢ１において、近傍半径ｒと、分類すべきクラスの個数ｋを設定する。近傍半径は、近傍半径内にあるすべての部分空間を同じクラスに属するとみなし、大雑把にクラス分類するための基準値として用いられる。例えばカラー画像の場合、同じような色のピクセル同士はＲＧＢ値も近く、したがって同じクラスに分類するのが自然であると考えられるからである。
【００７０】
後に説明するように、近傍半径が大きすぎる場合は所望のクラス数に到達する前に分類が完了してしまう場合があるが、近傍半径を縮小して再度クラス分類を行うため、最終的には必要な数のクラスに分類することができる。したがってこの近傍半径の初期値は十分大きな値に、例えば５０などに設定する。分類すべきクラスの個数ｋは、所与の値をそのまま用いる。
【００７１】
次に、ステップＢ２において、各部分空間ごとに、各部分空間の代表値Ｘｊを設定する。各部分空間の代表値としては、例えば部分空間の中央値などを用いることができる。以下では、これら代表値間の距離を部分空間の距離とする。
【００７２】
次に分類対象となる部分空間の集合を設定する（ステップＢ３）。以下ではこの集合をΩと書く。Ωの初期値は、データを含む全部分空間全体から成る集合とする。また分類済みクラスの数ｉを1とし、分類済みクラスＣｉの初期値を空集合とする。
【００７３】
次に、Ωに属する部分空間の中から、粗視化経験確率の最も高い部分空間を選択する（ステップＢ４）。この部分空間をＡＳと書く。次に、Ωに属する部分空間とＡＳとの距離を調べ、近傍半径ｒ以内にある部分空間を全て選択し、これを近傍集合とする（ステップＢ５）。以下では、近傍集合をＢＳと書く。
【００７４】
次に、すでに分類が完了しているクラスＣｉに含まれる部分空間と、近傍集合ＢＳに含まれる部分空間との最短距離を求め、近傍半径ｒよりも大きいかどうかを調べる（ステップＢ６）。もしこの最短距離がｒよりも大きいならば、近傍集合ＢＳは、すでに分類が完了したクラスと十分異なる属性を持ち、かつ高い確率で出現しているから、新しいクラスと見なしてよい。したがって、近傍集合ＢＳをそのまま新規クラスとして採用する。
【００７５】
ＢＳは分類が完了したので、分類対象集合Ωから削除する。図５においてはこの削除を、記号を用いて「Ω←Ω＼ＢＳ」と表している。Ωを更新したらステップＢ４に戻る（ステップＢ７）。
【００７６】
もし、すでに分類が完了しているクラスＣｉに含まれる部分空間と、近傍集合ＢＳに含まれる部分空間との最短距離が、近傍半径ｒよりも小さければ、近傍集合ＢＳはＣｉと属性が近いと考えられるため、ＢＳをＣｉに統合する。ＢＳは分類が完了したので、分類対象集合Ωから削除する（ステップＢ８）。
【００７７】
次に、Ωが空集合であるかどうかを調べ（ステップＢ９）、もし、空集合でなければ、ステップＢ４へ行き、Ωが空集合であれば、分類が完了したクラスの数がｋ個以上あるかどうかを調べ（ステップＢ１０）、ｋ個以下なら近傍半径ｒに１以下の定数をかけて半径を縮小する。この定数としては例えば0.9などの値を用いることができる。そののちステップＢ３以下を繰り返す。
【００７８】
もし、Ωが空集合であり、かつ分類済みのクラスの数が所与の数よりも多ければ、所望の数のクラス分類が完了しているので、各クラス内で前記クラスパラメータを算出し、これをクラスパラメータの初期値とし、また各クラス内に含まれる部分空間の数の比率を前記クラス混合比の初期値とする（ステップＢ１２）。
【００７９】
図７は、図５に示す粗視化経験確率分布に基づいてパラメータの初期値を決定する初期値設定の手順を理解しやすくするための簡単な例を表した図である。図７ではデータ空間が一次元であるとし、全部で１０個の部分空間に分けられていると仮定する。以下、図７を用いて初期値設定の手順について説明する。
【００８０】
図７（ａ）において、横軸は部分空間の番号を表し、縦軸は粗視化経験確率分布を表す。以下では処理の流れを直観的に説明することに主眼を置くため、粗視化経験確率値や部分空間の代表値、近傍半径については具体的な数値は用いず、図示することにする。
【００８１】
まず、ステップＢ１において、例えばクラス数を２とし、近傍半径をｒとする。ステップＢ２において、各部分空間の代表値を設定する。ステップＢ３において、分類対象集合Ωの初期値は、データを含んだ全部分空間とするので、
Ω＝｛Ａ３，Ａ４，Ａ５，Ａ６，Ａ７，Ａ８｝
となる。Ａ１，Ａ２およびＡ９，Ａ１０は、粗視化確率が０、すなわちこれらの部分空間に含まれるデータは観測されていないので、分類対象集合には含めない。
【００８２】
ステップＢ４で、分類対象集合に含まれる部分空間の中で、最も高い粗視化経験確率を持つＡ７を選び、これをＡＳとする（図７（ａ））。ステップＢ５で、Ａ７から近傍半径ｒ内にある部分空間を選択し、これをＢＳとする。図７（ａ）に示された近傍半径内にある部分空間は、Ａ５，Ａ６，Ａ７，Ａ８であるから、
ＢＳ＝｛Ａ５，Ａ６，Ａ７，Ａ８｝
となる（図７（ｂ））。
【００８３】
ステップＢ６において、分類されているクラスはまだないので、このＢＳをそのまま最初のクラスＣ１として採用し、分類対象集合からＢＳを取り除き、ステップＢ４に戻る。図７（ｂ）において、粗視化経験確率の高さを示す棒グラフが白抜きで表されているのは、部分空間が分類対象集合から除かれたことを表す。
【００８４】
ステップＢ４において、残りの分類対象集合の中で最も高い粗視化経験確率を持つものは、Ａ４であるから、これを新たにＡＳとする（図７（ｃ））。ステップＢ５で、Ａ４から近傍半径ｒ内にある部分空間を選択して、これをＢＳとする。ここでは、
ＢＳ＝｛Ａ３，Ａ４｝
となる。
【００８５】
ステップＢ６で、分類済みのクラス、すなわち
Ｃ１＝｛Ａ５，Ａ６，Ａ７，Ａ８｝
を調べると、Ａ４から近傍半径ｒ内にあるＡ５，Ａ６を含んでいる。したがって現在のＢＳを分類済みのクラスＣ１に統合する（図７（ｄ））。
【００８６】
これで分類対象集合Ωは空となり、全ての部分空間の分類が完了したことになるが、分類されたクラス数は１であり、所望のクラス数２に達していない（ステップＢ１０）。したがって近傍半径を縮小し（ステップ１１）、ステップＢ３以下を繰り返す。
【００８７】
以下、縮小された半径をｒ’とし（図７（ｅ））、上述の説明と同様の手続きを繰り返す。ただし、今回は近傍半径が縮小されているため、以下のような違いが生じる。即ち、Ａ７の近傍半径ｒ’内にある部分空間は、今回の場合
ＢＳ＝｛Ａ６，Ａ７，Ａ８｝
となる（図７（ｆ））。
【００８８】
このＢＳをそのまま最初のクラスＣ１として採用し、残りの分類対象集合の中から最も高い粗視化経験確率を持つＡ４を選ぶ（図７（ｇ））。Ａ４から近傍半径ｒ’内にある部分空間は、
ＢＳ＝｛Ａ３，Ａ４，Ａ５｝
となる。今回は、分類済みのクラスＣ１の中に、Ａ４から近傍半径ｒ’内にある部分空間は存在しないので、現在のＢＳを新規クラスＣ２として採用する（図７（ｈ））。これで全ての部分空間が所望の二つのクラスに分類されたことになる。
【００８９】
大雑把なクラス分類が完了すれば、分類されたクラス内で平均や分散を求め、以降で行う推定の初期パラメータとすることができる。初期パラメータを適切に設定することは、粗視化対数尤度の最大化の過程で局所最適解に陥ることを防止する上で有効である。
【００９０】
粗視化条件付き確率分布算出装置４では、このようにして求めたパラメータを初期値として粗視化条件付き確率分布を求めることが可能となる。粗視化条件付き確率分布算出装置４は、求めた粗視化条件付き確率分布をクラス帰属確率算出装置５に送る。クラス帰属確率算出装置５は、本発明第１記載の第５の工程に基づいて、式（１１）を用いてクラス帰属確率を算出する（ステップＡ６）。
【００９１】
このクラス帰属確率は、ｊ番目の部分空間に含まれるピクセルが、ｉ番目のクラスに属する確率を表す。したがって、画像各ピクセルに対してこのクラス帰属確率を算出し、確率の高いクラスに各ピクセルを分類することによって、領域抽出が容易に行うことが可能となる。クラス帰属確率算出装置５は、算出したクラス帰属確率をパラメータ更新装置６に送る。
【００９２】
パラメータ更新装置６では、本発明第１記載の第６の工程に基づいて、式（９）を最大化するようにパラメータを更新する。具体的には、以下のようにパラメータを更新する。
【００９３】
【数１６】
ここで、
【００９４】
【数１７】
は、ベクトルｕ，ｖのｉ成分とｊ成分の積ｕｉｖｊをｉｊ成分に持つマトリクスを表す。また、
【００９５】
【数１８】
は、式（１４）で定義されるｊ番目の部分空間Ａｊに含まれるデータの平均値である。
【００９６】
すでに述べたように、このようにパラメータを更新すれば、式（９）で与えられる粗視化された完全対数尤度の期待値は増加し、したがって式（１２）で与えられる粗視化平均対数尤度も増加する。このことは文献２に詳しい説明がある。パラメータ更新装置６は、パラメータを更新した後、更新したパラメータを評価関数算出装置７に送る。評価関数算出装置７は、本発明第１記載の第７の工程に従って、式（１２）を用いて粗視化対数尤度を算出する（ステップＡ８）。
【００９７】
また粗視化対数尤度のかわりに、本発明第４に記載のＡＩＣ
【００９８】
【数１９】
を用い、ＡＩＣが小さいほど良い推定結果であるとする評価関数を用いることも可能である。ここでｍは、全パラメータの個数を表す。ＡＩＣは粗視化対数尤度にマイナスを掛けた量に比例するため、パラメータ変更装置６で行う更新によってＡＩＣは減少する方向に変化する。また、パラメータに比例する項が付加されているため、同じ粗視化対数尤度ならばパラメータの少ないモデルを用いた推定結果の方が良いとする。この評価関数を用いることにより、データへの過剰な適合が抑制され、雑音に強い推定を行うことが可能となる。
【００９９】
また、本発明第５に記載するように、次式で定義されるＭＤＬ
【０１００】
【数２０】
を用いても同様な効果が得られる。
【０１０１】
さらに、ストラクチュラル・リスク
【０１０２】
【数２１】
を用いても同様な効果が得られる。ここでηは式（２１）が確率ηで成り立つことを表し、通常0.01などの値を用いる。c,a1,a2は確率分布の性質によって決まる定数であるが、通常
c＝１,a1＝１,a2＝１などの値が用いられる。ｈはＶＣ次元と呼ばれるものであり、パラメータの数に比例する量である。
【０１０３】
評価関数算出装置７は、本発明請求項１の第８の工程として記載されるように、評価関数の変化が所与の終了条件を満たしているかどうかを調べ、終了条件が満たされていれば領域抽出装置８に現在のパラメータを送り、終了条件が満たされていなければ粗視化条件付き確率分布算出装置４に現在のパラメータを送る（ステップＡ９）。終了条件としては、例えば現時点での評価関数値と前回算出した評価関数値の差を現時点での評価関数値で割り、その値の絶対値が0.0001以下であるかどうか、などの条件を用いることができる。
【０１０４】
領域抽出装置８は、本発明請求項１の第９の工程として記載されるように、評価関数算出装置７からパラメータを受け取り、パラメータの情報を用いて領域を抽出する（ステップＡ１０）。例えば、ｉ番目のクラスに属する領域を抽出するには、式（１１）で与えられる粗視化されたクラス帰属確率の値をｊ＝１からｊ＝Ｎまで調べ、その確率がある値（閾値）以上の部分空間をｉ番目のクラスに属する部分空間とする。次に、その部分空間に含まれるピクセルを調べ、これらのピクセルをｉ番目の領域であるとして抽出する。クラス帰属確率の閾値としては、例えば0.5を用いれば、所望の結果を得ることができる。
【０１０５】
領域抽出装置８は、本発明請求項４に記載されるように、自動的に閾値を設定することも可能である。そのためには以下のようにする。ｉ番目の領域を抽出するためには、まず推定されたクラス混合比ｗｉに全ピクセル数を掛け、各クラスに属するピクセル数の推定値を求める。この数をｎｉとする。
【０１０６】
次に、式（１１）の粗視化されたクラス帰属確率をｊ＝１からｊ＝Ｎまで調べ、値の高い部分空間から順に、その部分空間に含まれるピクセルを抽出し、抽出したピクセルがｎｉに到達するまで続ける。ｎｉ番目に抽出したピクセルの番号をｌとすると、式（１１）の粗視化されたクラス帰属確率の値が、ｉ番目の領域に属する確率の閾値となる。領域抽出装置８は、領域抽出が完了した後、抽出した領域のデータを出力装置９へ送る。
【０１０７】
上に説明した、データ空間粗視化装置２、粗視化経験分布算出装置３、粗視化条件付き確率分布４、クラス帰属確率算出装置５、パラメータ更新装置６、評価関数算出装置７、領域抽出装置８は、例えば、パーソナルコンピュータやワークステーション、あるいはスーパーコンピュータなどの計算機を用いて構築することができる。出力装置９は、領域抽出装置８から領域データを受け取り、出力する（ステップＡ１１）。出力装置９は、例えばディスプレイ装置やプリンタなどを用いて構成することができる。
【０１０８】
図３は本発明の第２の実施形態の画像領域抽出方法における処理手順を示すフローチャートである。また、図４は本実施形態の画像領域抽出装置の構成を示すブロック図である。本実施形態は、第１の実施形態に領域復元装置１０を追加した点を特徴とする。従って、以下の説明では、第１の実施の形態と重複する部分の説明は省略する。
【０１０９】
領域復元装置１０は、請求項３に記載されるように、粗視化確率分布を用いたパラメータ推定が完了した後、粗視化の解像度が元の解像度に等しいかを調べ（ステップＡ１２）、粗視化されていれば、データを元の解像度に戻す（ステップＡ１３）。粗視化されていなければ、元の解像度での推定が完了していることを意味するから、元の解像度で推定されたパラメータを領域抽出装置８に送る。
【０１１０】
データを元の解像度に戻すには、粗視化の解像度をデータ属性値の最小単位に設定し（例えば１）、第１の実施の形態に述べた方法と全く同一の方法を繰り返せばよい。この場合、第１の実施例で述べた方法よりも推定に時間を要するようになるが、より高い精度でパラメータを推定することができ、結果としてより精度の高い領域抽出が可能となる。また、粗視化確率分布を用いて推定したパラメータは、すでに最適なパラメータの近傍に推定されているため、最初から高い解像度で推定を行うよりも少ない回数のパラメータ更新で最適パラメータが推定でき、はるかに高速に領域抽出を行うことができる。
【０１１１】
次に、シミュレーションデータを用いて本発明の有効性を評価した結果について説明する。
【０１１２】
図８は、細胞の顕微鏡写真を模擬した画像であり、腺腔（画像中もっとも明るい領域）、細胞核（腺腔の周りを取り囲むように存在する輝度の低い小さな円形の組織）、およびそれ以外の組織から成る。書面ではモノクロ画像であるが、実際は腺腔が薄いピンク、細胞核が暗い青紫色、それ以外の領域が赤紫色のカラー画像となっている。
【０１１３】
ここでの目的は、画像から腺腔領域と細胞核領域を分離して抽出することである。抽出した腺腔の形状、細胞核の密度や大きさは、例えば細胞が癌化しているかどうかの判定に役立てることができる。しかしながら図８からわかるように、細胞核に対応する領域は背景と区別が付きにくく、通常の方法では細胞核領域だけを自動抽出することは困難である。
【０１１４】
この画像に対し、本発明に記載する手法を用いて、腺腔領域と細胞核領域を抽出することを試みた。ここでの推定は、粗視化解像度を16×16×16とし、クラスの数を３と仮定して行っている。クラスの数を３としたのは、細胞核に対応する領域（クラス１とする）、腺腔に対応する領域（クラス２とする）、それ以外の領域（クラス３とする）に分離できることを期待しているからである。
【０１１５】
このようにして得られた推定結果を図９（ａ）、（ｂ）、（ｃ）、（ｄ）に示す。図９（ａ）は、クラス１（細胞核領域）のクラス帰属確率を３次元的に表現したものであり、ｘ軸、ｙ軸が画像上のピクセル位置を示し、ｚ軸が帰属確率値を表している。図９（ｂ）はクラス１への帰属確率を２次元的に濃淡で表示したものであり、クラス１に属する確率が低いピクセルほど黒く、確率が高いピクセルほど白く表現されている。
【０１１６】
同様にして図９の（ｃ）、（ｄ）は、クラス２（腺腔領域）へのクラス帰属確率を、それぞれ３次元的、２次元的に表したものであり、帰属確率の高い部分が腺腔領域に属するピクセルとなっていることが分かる。（クラス３に対する結果は、今の場合抽出対象ではないので省略する。）これらの結果は、閾値を陽に設定することなく自動的に得られたものであり、各クラスへの帰属確率の分布を調べるだけで、領域抽出が可能となっている。
【０１１７】
同じ推定を、粗視化無し（粗視化解像度１）で行った結果を図１０の（ａ）、（ｂ）、（ｃ）、（ｄ）に示す。図１０（ａ）、（ｂ）は、クラス１へのクラス帰属確率を、それぞれ３次元的、２次元的に表したもの、また図１０（ｃ）、（ｄ）は、クラス２への帰属確率を、それぞれ３次元的、２次元的に表したものである。これらの図から、粗視化解像度１６での推定結果と粗視化無しでの推定結果がほとんど同じであることが分かる。
【０１１８】
推定されたパラメータ値を比較すると、例えばクラス１に属する領域の平均ＲＧＢ値は、粗視化無しで推定した場合ＲＧＢ＝（１４３.０６８，８６.６８９，１３５.７３７）、解像度１６で推定した場合ＲＧＢ＝（１４１.５２２，８６.１２８，１３５.５１３）となっており、粗視化の影響はほとんどない。
【０１１９】
しかしながら、粗視化のある無しでは推定時間に大きな差が生じる。この違いを調べるために、さまざまな粗視化解像度で推定したときの推定時間を比較したグラフを図１１に示す。このグラフより、粗視化解像度１６で推定した場合、粗視化無しで推定したときに比べて約１／７の時間で推定できていることがわかる。
【０１２０】
【発明の効果】
以上に説明したように本発明では、クラス帰属確率に基づいて領域抽出を行うことにより、閾値を陽に設定することなく、自動的に画像データから所望の領域を抽出することができる。また、粗視化確率分布を用いることにより、従来例よりも更に高速に領域抽出を行うことができる。
【図面の簡単な説明】
【図１】本発明の第１の実施の形態における処理手順を示すフローチャートである。
【図２】本発明の第１の実施の形態における、画像領域抽出装置の構成を示すブロック図である。
【図３】本発明の第２の実施の形態における処理手順を示すフローチャートである。
【図４】本発明の第２の実施の形態における、画像領域抽出装置の構成を示すブロック図である。
【図５】本発明の第１の実施の形態における、パラメータの初期化手順を示すフローチャートである。
【図６】粗視化データ空間と部分空間の一形態を示す図である。
【図７】 本発明請求項７記載の、パラメータの初期化方法の一例を表した図である。
【図８】本発明が適応できる画像の一例である。
【図９】本発明を用いた画像領域抽出を、粗視化解像度１６で行った結果を示しており、（ａ）クラス１（細胞核領域）に属するクラス帰属確率分布の３次元表示、（ｂ）クラス１（細胞核領域）に属するクラス帰属確率分布の２次元表示、（ｃ）クラス２（腺腔領域）に属するクラス帰属確率分布の３次元表示、（ｄ）クラス２（腺腔領域）に属するクラス帰属確率分布の２次元表示をそれぞれ示す。
【図１０】本発明を用いた画像領域抽出を、粗視化無し（粗視化解像度１）で行った結果を示しており、（ａ）クラス１（細胞核領域）に属するクラス帰属確率分布の３次元表示、（ｂ）クラス１（細胞核領域）に属するクラス帰属確率分布の２次元表示、（ｃ）クラス２（腺腔領域）に属するクラス帰属確率分布の３次元表示、（ｄ）クラス２（腺腔領域）に属するクラス帰属確率分布の２次元表示をそれぞれ示す。
【図１１】本発明を用いた画像領域抽出を、粗視化解像度１、２、４、８、１６で行ったときの推定時間を示すグラフである。
【符号の説明】
１ 入力装置
２ データ空間粗視化装置
３ 粗視化経験確率分布算出装置
４ 粗視化条件付き確率分布算出装置
５ クラス帰属確率算出装置
６ パラメータ更新装置
７ 評価関数算出装置
８ 領域抽出装置
９ 出力装置
１０ 領域復元装置
１１ 粗視化データ空間
１２ 部分空間
Ａ１〜Ａ１３ ステップ
Ｂ１〜Ｂ１２ ステップ
【発明の属する技術分野】
本発明は、領域抽出の対象像を含む画像データに対し、画像のピクセルが持つ属性に基づいて、対象となる物体の領域を抽出する技術に関し、特に、ＭＲＩ装置やＣＴスキャン装置などで撮影された画像データ、あるいは生物顕微鏡で観測された生体組織の断面図から、画像を特定の臓器や腫瘍、細胞の核、腺腔などの領域を抽出する画像領域抽出方法および装置に関する。
【０００２】
【従来の技術】
従来、画像から特定の物体の像を抽出するには、次のような方法が知られている（例えば、特開２００１−９２９８０号公報；以下文献１とする。）
文献１にでは、輪郭抽出の対象となる物体の像を含む画像データに対し、前記画像の各点が持つ属性に基づいて、画像各点が各領域に属する確率である領域帰属確率を用いて、前記各点が属する領域を分離し、前記領域間の境界を輪郭として抽出する構成としている。このような構成とすることにより、領域分類の閾値を陽（explicit）に設定することなく、自動的に各領域の輪郭を抽出することが可能となり、輪郭領域抽出の高速化を可能にしている。
【０００３】
【発明が解決しようとする課題】
上記文献１における輪郭領域抽出は、画像データから、画像上各点が領域の各々に属する領域帰属確率の期待値を算出し、領域帰属確率と領域パラメータから求められる混合確率分布に基づいて評価関数を算出し、領域帰属確率に基づいて各領域を分離し、分離された領域に基づいて輪郭を抽出しているため、評価関数を求める際、すべての画像点（ピクセル）に関して和を取る必要があり、最適なパラメータを求めるためにはこの評価関数の算出を何度も繰り返すことになる。
【０００４】
そのため、従来の方法においては，画像のサイズが大きくなればなるほど、領域抽出に要する時間が膨大になるという問題がある。
【０００５】
本発明は上述した問題点を鑑み、画像から、閾値を陽に設定することなく、更なる高速かつ高精度に対象領域を抽出する画像領域抽出方法および装置を提供することにある。
【０００６】
【課題を解決するための手段】
本発明第１記載の画像領域抽出方法は、
領域抽出の対象となる画像データに対し、画像上の各ピクセルを複数のクラスに分割し、所望の領域を抽出する画像領域抽出方法において、
前記画像の前記各ピクセルが取りうる属性値全体から構成されるデータ空間を、所与の解像度で部分空間に分解し、該各部分空間に属性値を取るピクセルの番号と、該ピクセルの属性値の平均値と、該ピクセルの個数とを保持し、粗視化データ空間を構成する第１の工程と、
前記各部分空間のピクセルの個数を前記画像に含まれる全ピクセルの個数で割り、前記粗視化データ空間上の粗視化経験確率分布を算出する第２の工程と、
前記各クラスの属性を規定するクラスパラメータ、前記クラスの数、および前記クラスの混合比率を初期化する第３の工程と、
前記クラスパラメータに基づいて、前記各部分空間の前記平均値における条件付き確率を算出し、該条件付き確率のデータ空間における分布である粗視化条件付き確率分布を算出する第４の工程と、
前記画像を構成する各ピクセルが各前記クラスに帰属する確率であるクラス帰属確率を、前記混合比率に前記粗視化条件付き確率分布を掛けることにより算出する第５の工程と、
粗視化対数尤度を評価関数として用い、該評価関数を増加させるように前記クラスパラメータ、および前記クラスの混合比率を更新する第６の工程と、
前記粗視化条件付き確率分布と、前記粗視化経験確率分布と、前記第６の工程で更新した前記クラスパラメータおよび前記クラスの混合比率とを用いて、前記評価関数を算出する第７の工程と、
前記評価関数が所与の終了条件を満たしているかどうかを調べる第８の工程と、
前記評価関数が前記所与の終了条件を満たした後、前記クラスパラメータおよび前記混合比率を保持し、前記クラス帰属確率に基づいて、前記ピクセルが属するクラスを決定し、前記所望の領域を抽出する第９の工程とを有し、
前記評価関数が前記所与の条件を満たすまで、前記第４、前記第５、前記第６、前記第７、および前記第８の工程を繰り返すことを特徴とする。
【０００８】
本発明第２に記載の画像領域抽出方法は、
前記第１に記載の画像領域抽出方法において、粗視化の解像度が元の解像度であるかを調べる第１０の工程と、前記部分空間の解像度を元の解像度へ復元する第１１の工程とを有し、
前記第３の工程において、前記第９の工程で保持されているクラスパラメータおよび前記クラスの混合比率を初期値として用い、所与の条件が満たされるまで、前記第４、前記第５、前記第６、前記第７、および前記第８の工程を繰り返すことを特徴とする。
【０００９】
本発明第３記載の画像領域抽出方法は、
前記第１及至第２のいずれかに記載の画像領域抽出方法における前記第９の工程において、推定された前記混合比率に前記画像を構成する全ピクセルの個数を掛けて前記各クラスに属するピクセルの個数を算出し、前記クラス帰属確率が高い順から前記ピクセルを選択し、前記各クラスに属するピクセルを決定することを特徴とする。
【００１０】
本発明第４記載の画像領域抽出方法は、
前記第１及至第３のいずれかに記載の画像領域抽出方法における前記第６の工程において、前記評価関数としてＡＩＣ（Akaike Information Criterion）を用い、前記第７の工程において評価関数を減少させるようパラメータを変更することを特徴とする。
【００１１】
本発明第５記載の画像領域抽出方法は、
前記第１及至第４のいずれかに記載の画像領域抽出方法における前記第６の工程において、前記評価関数としてＭＤＬ（Minimum Description Length）を用い、前記第７の工程において評価関数を減少させるようパラメータを変更することを特徴とする。
【００１２】
本発明第６記載の画像領域抽出方法は、
前記第１及至第４のいずれかに記載の画像領域抽出方法における前記第６の工程において、前記評価関数としてストラクチュラル・リスク（Structural Risk）を用い、前記第７の工程において該評価関数を減少させるようパラメータを変更することを特徴とする。
【００１３】
本発明第７記載の画像領域抽出方法は、
前記第１及至第６のいずれかに記載の画像領域抽出方法における前記第３の工程において、
前記各部分空間同士が近傍にあるかどうかを定義する近傍半径と、前記クラスの数を設定し（ステップＢ１）、
前記各部分空間ごとに、各部分空間の代表値を設定し（ステップＢ２）、
分類対象部分空間の集合を設定し（ステップＢ３）、
前記分類対象部分空間の集合の中から、前記粗視化経験確率の最も高い部分空間を選択し（ステップＢ４）、
前記粗視化経験確率の最も高い部分空間の代表値との距離が、前記近傍半径以内である代表値を持つ部分空間全てを、近傍集合として選択し（ステップＢ５）、
すでに分類が完了したクラスに含まれる部分空間の代表値と、前記近傍集合に含まれる部分空間の代表値との間の最短距離が、前記近傍半径よりも大きいかどうかを調べ（ステップＢ６）、
もし、すでに分類が完了したクラスに含まれる部分空間の代表値と、前記近傍集合に含まれる部分空間の代表値との間の最短距離が、前記近傍半径よりも大きければ、前記近傍集合を新規クラスとし、前記分類対象部分空間から前記近傍集合を削除し、ステップＢ４以降を繰り返し（ステップＢ７）、
そうでなければ、前記近傍集合をすでに分類済みのクラスに追加し、前記分類対象部分空間から、前記近傍集合を削除し（ステップＢ８）、
前記分類対象部分空間が空集合であるかどうかを調べ（ステップＢ９）、
もし、前記分類対象部分空間が空集合でなければ、ステップＢ４以下を繰り返し、
前記分類対象部分空間が空集合であれば、分類が完了したクラスの数が所与の数以上あるかどうかを調べ（ステップＢ１０）、
もし分類が完了したクラスの数が所与の数よりも少なければ、前記近傍半径を縮小し（ステップＢ１１）、ステップＢ３以下を繰り返し、
もし、前記分類対象部分空間が空集合であり、かつ分類済みのクラスの数が所与の数よりも多ければ、各クラス内で前記クラスパラメータを算出し、これをクラスパラメータの初期値とし、また各クラス内に含まれる部分空間の数の比率を前記混合比率の初期値とする（ステップＢ１２）、ことを特徴とする。
【００１４】
本発明第８記載の画像領域抽出装置は、
領域抽出の対象となる像を含む画像データに対し、画像上の各ピクセルを複数のクラスに分類し、所望の領域を抽出する画像領域抽出装置において、
画像データを読み込む入力装置と、
前記画像の前記各ピクセルが取りうる属性値全体から構成されるデータ空間を、所与の解像度で部分空間に分解し、該各部分空間に属性値を取るピクセルの番号と、該ピクセルの属性値の平均と、該ピクセルの個数を保持し、粗視化データ空間を構成する領域粗視化装置と、
前記各部分空間のピクセルの個数を前記画像に含まれる全ピクセルの個数で割り、前記粗視化データ空間上の粗視化経験分布を算出する粗視化経験確率分布算出装置と、
前記各クラスの属性を規定するクラスパラメータ、前記クラスの数、および前記クラスの混合比率を初期化し、前記各クラスの属性を規定するクラスパラメータから、前記クラスが指定されたもとでの条件付き確率分布を算出し、前記クラスが指定されたもとでの条件付き確率分布を、前記各部分空間内で平均化することによって粗視化条件付き確率分布を算出する粗視化条件付き確率分布装置と、
前記画像を構成する各ピクセルが前記各クラスに帰属する確率であるクラス帰属確率を、前記粗視化条件付き確率分布から算出するクラス帰属確率算出装置と、
評価関数として粗視化対数尤度を用いるときには、該評価関数を増加させるように、また評価関数としてＡＩＣ、ＭＤＬ、ストラクチュラル・リスクのいずれかを用いる場合は該評価関数を減少させるように前記パラメータを更新するパラメータ更新装置と、
前記粗視化対数尤度、前記ＡＩＣ、前記ＭＤＬ、前記ストラクチュラル・リスクのいずれかを用いて前記評価関数を算出する評価関数算出装置と、
前記評価関数が所与の終了条件を満たしているかどうかを調べ、前記評価関数が前記所与の終了条件を満たした後、前記パラメータを保持し、前記クラス帰属確率に基づいて、前記各ピクセルが属するクラスを決定し、領域を抽出する領域抽出装置と、
抽出した領域を出力する出力装置と、
を有することを特徴とする。
【００１５】
本発明第９記載の画像領域抽出装置は、前記第８記載の画像領域抽出装置において、前記評価関数が前記所与の終了条件を満たしていることが確認されたのち、粗視化の解像度が元の解像度であるかを調べ、前記データ空間の解像度を元の解像度へ復元する解像度復元装置を有することを特徴とする。
【００１６】
また、本発明は、本発明における画像領域抽出処理をコンピュータに実行させるための画像領域抽出プログラムをその特徴とする。
【００１７】
本発明では文献１記載の発明と同様に、推定したクラス帰属確率に基づいて領域を抽出している。この特徴により、閾値を自動的に設定できるという利点がある。しかしながら、このクラス帰属確率を推定するためには、パラメータの更新を何度も（通常５０回程度）繰り返し行う必要があり、計算に時間を要するという問題がある。
【００１８】
そこで、本発明では画像データの発生源のモデルとして粗視化確率分布を用い、それに伴ってデータ空間も粗視化する。これにより、通常は条件付き確率分布やクラス帰属確率を数十万個のデータについて計算する必要があるところを、数千個の部分空間についての計算に置き換えることができる。この特徴により、領域抽出に要する時間を大幅に短縮することが可能となる。
【００１９】
推定されるクラスパラメータは、粗視化によって精度が劣ることになるが、本発明ではクラスパラメータを直接使うのではなく、パラメータから計算したクラス帰属確率に基づいて領域抽出を行うため、後に示すように、粗視化による誤差の影響をうけることなく領域抽出を行うことが可能となっている。
【００２０】
以下、上記本発明の作用の有効性について更に詳細に説明する。
【００２１】
本発明の第１に記載の画像領域抽出方法においては、画像を構成する各ピクセルの持つ属性値を確率変数と見なし、推定したピクセル値の確率分布に基づいて所望の領域を抽出する。属性値としては、例えばモノクロ画像ならば輝度の値、カラー画像ならば赤（Ｒ）、緑（Ｇ）、青（Ｂ）の色要素の強度などを用いることができる。
【００２２】
所望の領域を抽出するためには、各ピクセルの属性値に基づいて、類似の属性をもった複数のグループに各ピクセルを分類する必要がある。本明細書においては、類似の属性を持ったピクセルの集合をクラスと呼ぶ。各クラスは、例えばそのクラスに属するピクセルの属性値の平均値や分散などによって特徴づけられるものとする。以下、これらクラスの特徴をそのクラスの「クラスパラメータ」と呼び、φi（i=1,,,k）と書くこととする。ｋはクラスの数である。このような定式化のもとで、ｊ番目のピクセルが、ｘjという属性値を取る確率は、次の混合分布で表現できる。
【００２３】
【数１】
【００２４】
【数２】
【００２５】
ｘjは、例えばモノクロ画像ならば、0から255などの整数値をとり、カラー画像ならばＲＧＢの色要素の値を成分とした３次元ベクトル（xj1,xj2,xj3）で表される。ここで各ｘjl（l=1,2,3）は、例えば0から255の整数値を取る。
【００２６】
混合比ｗiは、異なるクラスに属する領域の面積比を表す。例えば、平均輝度200、輝度の標準偏差が20という特徴で表される明るい画像領域（クラス１とする）と、平均輝度50、輝度の標準偏差が10という特徴で表される暗い領域（クラス２とする）の、二つの領域から構成されるモノクロ画像があるとする。また、その画像の面積の７割を明るい領域が、３割を暗い領域が占めているとする。この場合、クラスの数はｋ=2、クラスパラメータはφ1=(200,20),φ2=(50,10)であり、この画像の混合分布は、
p(xj)=0.7(xj|200,20)+0.3(xj|50,10) ……(２）
などと表現できる。
【００２７】
以下、クラス混合比ｗiとクラスパラメータφiをまとめてθiと書くことにする。以下において単に「パラメータ」というときは、このθiを意味することとする。
【００２８】
本発明第１記載の画像領域抽出方法では、以下のように定義される平均対数尤度
【００２９】
【数３】
【００３０】
しかしながら、平均対数尤度Ｌを最大にするパラメータを推定することは、一般に困難である。そこで、この平均対数尤度Ｌの代わりに次の量
【００３１】
【数４】
【００３２】
【数５】
【００３３】
本発明において、推定されたパラメータから実際に領域を抽出するには、以下のように行う。まず適当な初期パラメータから開始し、（５）式で与えられるクラス帰属確率を求める。次にＱを増加させるようにパラメータｗ，φを更新し、改めてＱを算出する。最終的にＱが増加しなくなるまでこの手順を繰り返す。
【００３４】
ｋ個のクラスのうち、例えばｉ番目のクラスに属する点を抽出するには、ｉ番目のクラスへの帰属確率の値を調べ、その確率値が一定値以上のピクセルをそのクラスに属するものとして分類する。すべてのピクセルを分類し終わった時点で、ｋ個のクラスの中から所望の属性を持ったクラスを選択し、そのクラスに属するピクセルを画像から抽出すれば、所望の領域を自動的に抽出することができる。
【００３５】
本発明では、このＱの最大化を高速に行うために、粗視化確率分布を導入する。粗視化確率分布は、データが取りうる値全体からなる空間（以下、データ空間と呼ぶ）を互いに交わらないＮ個の部分空間に分解し、各部分空間に確率を割り当てることによって構成される。
【００３６】
具体的には、ｊ番目の部分空間における粗視化された条件付き確率分布を
【００３７】
【数６】
【００３８】
【数７】
【００３９】
【数８】
このように定義された粗視化条件付き確率分布を用いると、Ｑは、
【００４０】
【数９】
【００４１】
【数１０】
【００４２】
【数１１】
【００４３】
【数１２】
【００４４】
元のＱと比較すると、式（４）で与えられるＱが全データに関して和をとっているのに対し、式（９）で与えられる粗視化された完全対数尤度は、部分空間に関してのみ和を取っている。後に示すように、この改善によって本発明では計算量が大幅に軽減できる。例えば、512×512ピクセルの画像の場合、２６万個以上のデータについて和を取る必要があるが、本発明の粗視化分布を用いた場合、部分空間に関する和を１０００個程度にまで軽減でき、高速な推定が可能となる。
【００４５】
上記の方法において各部分空間における粗視化された確率値は、その部分空間に含まれるピクセルの属性値の平均値における確率値で近似する。
【００４６】
【数１３】
【００４７】
【数１４】
【００４８】
本発明第２に記載の画像領域抽出方法では、粗視化確率分布を用いて推定したパラメータを初期値とし、元の解像度で再度推定を行う。この場合、粗視化確率分布を用いてすでにほぼ最適なパラメータが得られているため、最初から元の解像度で推定を行う場合に比較して、パラメータの逐次更新の回数がはるかに少なくてすむため、高精度な推定を高速に行うことが可能となる。
【００４９】
本発明第３に記載の画像領域抽出方法では、領域を抽出する際、推定された混合比ｗiに全ピクセル数ｎを掛け、ｉ番目の領域に属するピクセル数を推定する。そして領域帰属確率の高い順に上位ｎi個をこの領域に属するピクセルとして抽出する。この方法により、どの値までの確率値をその領域に属しているものと見なすか、という閾値を自動的に決定することができる。
【００５０】
本発明第４、第５および第６に記載の画像領域抽出方法においては、評価関数として、それぞれＡＩＣ、ＭＤＬ、ストラクチュラル・リスクを用い、その値がもっとも小さい結果を与えるモデルを選択する。最尤推定においては一般にパラメータの数が多いほど、評価値としてはよい値を与えるため、本来のパラメータ数が決定できないという問題がある。
【００５１】
これに対し、本発明第４、第５および第６に記載の評価関数は、過剰なパラメータ数を用いたとき、評価関数の値が逆に増加するという特性を持っている。これにより、最適なパラメータ数を推定することが可能となる。領域抽出へ適用する場合は、画像が何種類の領域から形成されているか、その適切な領域数を推定することができる。
【００５２】
【発明の実施の形態】
次に、本発明の好ましい実施の形態について、図面を参照しながら詳細に説明する。なお、以下の説明において用いる記号は、特に言及しない限りこれまで用いてきた記号の説明に従うものである。また、以下の実施例では、染色された細胞のカラー顕微鏡写真から、細胞の各組織（核、腺腔）を抽出する場合を例に取って説明を行うが、画像として記録されている像であるならば、以下に説明する方法と同様な方法で実行できる。
【００５３】
図１は本発明の第１の実施形態の画像領域抽出方法における処理手順を示すフローチャートである。また、図２は本実施形態の画像領域抽出装置の構成を示すブロック図である。
【００５４】
本発明の画像領域抽出装置は、画像データを読み込む入力装置１と、データ空間粗視化装置２と、粗視化経験確率分布算出装置３と、粗視化条件付き確率分布算出装置４と、クラス帰属確率算出装置５と、パラメータ更新装置６と、評価関数算出装置７と、領域抽出装置８と、出力装置９から構成されている。
【００５５】
以下、本発明の画像領域抽出装置の処理について、図１および図２を用いて説明する。
【００５６】
入力装置１は、生物顕微鏡などで観測された細胞のカラー画像を入力する（ステップＡ１）。この入力装置は、例えば、画像スキャナやデジタルカメラなどを用いて構成することができる。あるいは、生物顕微鏡をコンピュータに接続し、ネットワークを通して画像を直接入力することも可能である。入力装置１は、読み込んだデータをデータ空間粗視化装置２に送る。
【００５７】
データ空間粗視化装置２は、本発明第１記載の第１の工程に基づいて、粗視化データ空間を構成する（ステップＡ２）。ここで、データ空間とは、画像の各ピクセルが取りうる属性値全体の集合を意味する。例えば標準的グレースケールの画像において、各ピクセルの属性は輝度によって表すことができ、輝度は通常0から255の整数値で表現される１次元空間である。また、標準的なカラー画像の場合ならば、Ｒ、Ｇ、Ｂ各色要素に対し、通常0から255の整数値で表現される３次元空間となる。
【００５８】
部分空間を構成するには、例えば解像度を８とすると、ＲＧＢの各値を８ずつに区切り、８×８×８の立方体を１つの部分空間として構成することができる。解像度は３次元の軸で同一で無くともよい。例えばＲ軸をｈ1ずつ、Ｇ軸をｈ2ずつ、Ｂ軸をｈ3ずつに区切って粗視化してもよく、各部分空間が重なりを持たず、かつ全データ空間をもれなく覆うことができればよい。以下において、解像度ｈで粗視化すると表現する場合は、ｈ×ｈ×ｈでＲＧＢの各値を区切ることを意味する。
【００５９】
データ空間粗視化装置２は、また、各部分空間に値を取るピクセルの集合と、それらピクセルの属性値の平均と、ピクセルの個数を保持する。
【００６０】
図６に解像度ｈ1×ｈ2×ｈ3で粗視化したときの粗視化データ空間１１、および部分空間１２を示す。例えばｊ番目のピクセルの属性値ｘjが、Ｒ=128、Ｇ=255、Ｂ=18なる値の時、このピクセルは、（16,31,2）なる指標で指定される部分空間に値を取る。このことを以下では、ｘjがこの部分空間に含まれる、と表現する。
【００６１】
データ空間粗視化装置２は、粗視化データ空間を構成した後、各部分空間に含まれるピクセルの番号と、ピクセル数を粗視化経験確率分布算出装置３に送る。粗視化経験確率分布算出装置３では、本発明第１記載の第２の工程に基づいて、粗視化経験確率分布を算出する（ステップＡ３）。ここで粗視化経験確率分布とは、粗視化データ空間の各部分空間に含まれるピクセルの個数を、全ピクセル数で割った値から算出された確率分布を表す。
【００６２】
粗視化経験確率分布算出装置３は、粗視化経験確率分布を粗視化条件付き確率分布算出装置４に送る。粗視化条件付き確率分布算出装置４では、本発明第１記載の第３の工程に基づいて、クラスパラメータを初期化する（ステップＡ４）。また、粗視化条件付き確率分布算出装置４では、本発明第１記載の第４の工程に基づいて、粗視化条件付き確率分布を算出する（ステップＡ５）。
【００６３】
粗視化条件付き確率分布を具体的に算出するには、以下のようにする。ここでは、ｊ番目のピクセル値がｉ番目のクラスから発生しているとした条件のもとでの条件付き確率が、以下のような多次元正規分布
【００６４】
【数１５】
【００６５】
粗視化条件付き確率分布算出装置４では、この条件付き確率を式（６）で与えられる式で計算する。このとき、各部分空間の測度ｍ(Ａj)は、各部分空間の体積となる。例えば、解像度８で一様に粗視化した場合、8×8×8=512となる。
【００６６】
粗視化条件付き確率分布算出装置４では、本発明第１記載の方法のように、粗視化条件付き確率分布を式（１３）によって近似することも可能である。この方法により、パラメータを更新する度に式（６）で与えられる演算をする必要がなくなり、計算量を大幅に削減できる効果がある。
【００６７】
また、粗視化条件付き確率分布算出装置４では、本発明第７記載の方法のように、粗視化経験分布に基づいてパラメータの初期値を決定することも可能である。本発明第７記載の方法では、各部分空間を大雑把にクラス分類し、分類の結果得られた各クラス内で平均値や分散値を求め、これらの値をパラメータ推定の初期値とする。
【００６８】
図５は、粗視化経験確率分布に基づいて、パラメータの初期値を決定する処理手順を示したフローチャートである。以下、図５を参照しながらパラメータの初期値設定法について説明する。
【００６９】
まず、ステップＢ１において、近傍半径ｒと、分類すべきクラスの個数ｋを設定する。近傍半径は、近傍半径内にあるすべての部分空間を同じクラスに属するとみなし、大雑把にクラス分類するための基準値として用いられる。例えばカラー画像の場合、同じような色のピクセル同士はＲＧＢ値も近く、したがって同じクラスに分類するのが自然であると考えられるからである。
【００７０】
後に説明するように、近傍半径が大きすぎる場合は所望のクラス数に到達する前に分類が完了してしまう場合があるが、近傍半径を縮小して再度クラス分類を行うため、最終的には必要な数のクラスに分類することができる。したがってこの近傍半径の初期値は十分大きな値に、例えば５０などに設定する。分類すべきクラスの個数ｋは、所与の値をそのまま用いる。
【００７１】
次に、ステップＢ２において、各部分空間ごとに、各部分空間の代表値Ｘｊを設定する。各部分空間の代表値としては、例えば部分空間の中央値などを用いることができる。以下では、これら代表値間の距離を部分空間の距離とする。
【００７２】
次に分類対象となる部分空間の集合を設定する（ステップＢ３）。以下ではこの集合をΩと書く。Ωの初期値は、データを含む全部分空間全体から成る集合とする。また分類済みクラスの数ｉを1とし、分類済みクラスＣｉの初期値を空集合とする。
【００７３】
次に、Ωに属する部分空間の中から、粗視化経験確率の最も高い部分空間を選択する（ステップＢ４）。この部分空間をＡＳと書く。次に、Ωに属する部分空間とＡＳとの距離を調べ、近傍半径ｒ以内にある部分空間を全て選択し、これを近傍集合とする（ステップＢ５）。以下では、近傍集合をＢＳと書く。
【００７４】
次に、すでに分類が完了しているクラスＣｉに含まれる部分空間と、近傍集合ＢＳに含まれる部分空間との最短距離を求め、近傍半径ｒよりも大きいかどうかを調べる（ステップＢ６）。もしこの最短距離がｒよりも大きいならば、近傍集合ＢＳは、すでに分類が完了したクラスと十分異なる属性を持ち、かつ高い確率で出現しているから、新しいクラスと見なしてよい。したがって、近傍集合ＢＳをそのまま新規クラスとして採用する。
【００７５】
ＢＳは分類が完了したので、分類対象集合Ωから削除する。図５においてはこの削除を、記号を用いて「Ω←Ω＼ＢＳ」と表している。Ωを更新したらステップＢ４に戻る（ステップＢ７）。
【００７６】
もし、すでに分類が完了しているクラスＣｉに含まれる部分空間と、近傍集合ＢＳに含まれる部分空間との最短距離が、近傍半径ｒよりも小さければ、近傍集合ＢＳはＣｉと属性が近いと考えられるため、ＢＳをＣｉに統合する。ＢＳは分類が完了したので、分類対象集合Ωから削除する（ステップＢ８）。
【００７７】
次に、Ωが空集合であるかどうかを調べ（ステップＢ９）、もし、空集合でなければ、ステップＢ４へ行き、Ωが空集合であれば、分類が完了したクラスの数がｋ個以上あるかどうかを調べ（ステップＢ１０）、ｋ個以下なら近傍半径ｒに１以下の定数をかけて半径を縮小する。この定数としては例えば0.9などの値を用いることができる。そののちステップＢ３以下を繰り返す。
【００７８】
もし、Ωが空集合であり、かつ分類済みのクラスの数が所与の数よりも多ければ、所望の数のクラス分類が完了しているので、各クラス内で前記クラスパラメータを算出し、これをクラスパラメータの初期値とし、また各クラス内に含まれる部分空間の数の比率を前記クラス混合比の初期値とする（ステップＢ１２）。
【００７９】
図７は、図５に示す粗視化経験確率分布に基づいてパラメータの初期値を決定する初期値設定の手順を理解しやすくするための簡単な例を表した図である。図７ではデータ空間が一次元であるとし、全部で１０個の部分空間に分けられていると仮定する。以下、図７を用いて初期値設定の手順について説明する。
【００８０】
図７（ａ）において、横軸は部分空間の番号を表し、縦軸は粗視化経験確率分布を表す。以下では処理の流れを直観的に説明することに主眼を置くため、粗視化経験確率値や部分空間の代表値、近傍半径については具体的な数値は用いず、図示することにする。
【００８１】
まず、ステップＢ１において、例えばクラス数を２とし、近傍半径をｒとする。ステップＢ２において、各部分空間の代表値を設定する。ステップＢ３において、分類対象集合Ωの初期値は、データを含んだ全部分空間とするので、
Ω＝｛Ａ３，Ａ４，Ａ５，Ａ６，Ａ７，Ａ８｝
となる。Ａ１，Ａ２およびＡ９，Ａ１０は、粗視化確率が０、すなわちこれらの部分空間に含まれるデータは観測されていないので、分類対象集合には含めない。
【００８２】
ステップＢ４で、分類対象集合に含まれる部分空間の中で、最も高い粗視化経験確率を持つＡ７を選び、これをＡＳとする（図７（ａ））。ステップＢ５で、Ａ７から近傍半径ｒ内にある部分空間を選択し、これをＢＳとする。図７（ａ）に示された近傍半径内にある部分空間は、Ａ５，Ａ６，Ａ７，Ａ８であるから、
ＢＳ＝｛Ａ５，Ａ６，Ａ７，Ａ８｝
となる（図７（ｂ））。
【００８３】
ステップＢ６において、分類されているクラスはまだないので、このＢＳをそのまま最初のクラスＣ１として採用し、分類対象集合からＢＳを取り除き、ステップＢ４に戻る。図７（ｂ）において、粗視化経験確率の高さを示す棒グラフが白抜きで表されているのは、部分空間が分類対象集合から除かれたことを表す。
【００８４】
ステップＢ４において、残りの分類対象集合の中で最も高い粗視化経験確率を持つものは、Ａ４であるから、これを新たにＡＳとする（図７（ｃ））。ステップＢ５で、Ａ４から近傍半径ｒ内にある部分空間を選択して、これをＢＳとする。ここでは、
ＢＳ＝｛Ａ３，Ａ４｝
となる。
【００８５】
ステップＢ６で、分類済みのクラス、すなわち
Ｃ１＝｛Ａ５，Ａ６，Ａ７，Ａ８｝
を調べると、Ａ４から近傍半径ｒ内にあるＡ５，Ａ６を含んでいる。したがって現在のＢＳを分類済みのクラスＣ１に統合する（図７（ｄ））。
【００８６】
これで分類対象集合Ωは空となり、全ての部分空間の分類が完了したことになるが、分類されたクラス数は１であり、所望のクラス数２に達していない（ステップＢ１０）。したがって近傍半径を縮小し（ステップ１１）、ステップＢ３以下を繰り返す。
【００８７】
以下、縮小された半径をｒ’とし（図７（ｅ））、上述の説明と同様の手続きを繰り返す。ただし、今回は近傍半径が縮小されているため、以下のような違いが生じる。即ち、Ａ７の近傍半径ｒ’内にある部分空間は、今回の場合
ＢＳ＝｛Ａ６，Ａ７，Ａ８｝
となる（図７（ｆ））。
【００８８】
このＢＳをそのまま最初のクラスＣ１として採用し、残りの分類対象集合の中から最も高い粗視化経験確率を持つＡ４を選ぶ（図７（ｇ））。Ａ４から近傍半径ｒ’内にある部分空間は、
ＢＳ＝｛Ａ３，Ａ４，Ａ５｝
となる。今回は、分類済みのクラスＣ１の中に、Ａ４から近傍半径ｒ’内にある部分空間は存在しないので、現在のＢＳを新規クラスＣ２として採用する（図７（ｈ））。これで全ての部分空間が所望の二つのクラスに分類されたことになる。
【００８９】
大雑把なクラス分類が完了すれば、分類されたクラス内で平均や分散を求め、以降で行う推定の初期パラメータとすることができる。初期パラメータを適切に設定することは、粗視化対数尤度の最大化の過程で局所最適解に陥ることを防止する上で有効である。
【００９０】
粗視化条件付き確率分布算出装置４では、このようにして求めたパラメータを初期値として粗視化条件付き確率分布を求めることが可能となる。粗視化条件付き確率分布算出装置４は、求めた粗視化条件付き確率分布をクラス帰属確率算出装置５に送る。クラス帰属確率算出装置５は、本発明第１記載の第５の工程に基づいて、式（１１）を用いてクラス帰属確率を算出する（ステップＡ６）。
【００９１】
このクラス帰属確率は、ｊ番目の部分空間に含まれるピクセルが、ｉ番目のクラスに属する確率を表す。したがって、画像各ピクセルに対してこのクラス帰属確率を算出し、確率の高いクラスに各ピクセルを分類することによって、領域抽出が容易に行うことが可能となる。クラス帰属確率算出装置５は、算出したクラス帰属確率をパラメータ更新装置６に送る。
【００９２】
パラメータ更新装置６では、本発明第１記載の第６の工程に基づいて、式（９）を最大化するようにパラメータを更新する。具体的には、以下のようにパラメータを更新する。
【００９３】
【数１６】
【００９４】
【数１７】
【００９５】
【数１８】
【００９６】
すでに述べたように、このようにパラメータを更新すれば、式（９）で与えられる粗視化された完全対数尤度の期待値は増加し、したがって式（１２）で与えられる粗視化平均対数尤度も増加する。このことは文献２に詳しい説明がある。パラメータ更新装置６は、パラメータを更新した後、更新したパラメータを評価関数算出装置７に送る。評価関数算出装置７は、本発明第１記載の第７の工程に従って、式（１２）を用いて粗視化対数尤度を算出する（ステップＡ８）。
【００９７】
また粗視化対数尤度のかわりに、本発明第４に記載のＡＩＣ
【００９８】
【数１９】
【００９９】
また、本発明第５に記載するように、次式で定義されるＭＤＬ
【０１００】
【数２０】
【０１０１】
さらに、ストラクチュラル・リスク
【０１０２】
【数２１】
c＝１,a1＝１,a2＝１などの値が用いられる。ｈはＶＣ次元と呼ばれるものであり、パラメータの数に比例する量である。
【０１０３】
評価関数算出装置７は、本発明請求項１の第８の工程として記載されるように、評価関数の変化が所与の終了条件を満たしているかどうかを調べ、終了条件が満たされていれば領域抽出装置８に現在のパラメータを送り、終了条件が満たされていなければ粗視化条件付き確率分布算出装置４に現在のパラメータを送る（ステップＡ９）。終了条件としては、例えば現時点での評価関数値と前回算出した評価関数値の差を現時点での評価関数値で割り、その値の絶対値が0.0001以下であるかどうか、などの条件を用いることができる。
【０１０４】
領域抽出装置８は、本発明請求項１の第９の工程として記載されるように、評価関数算出装置７からパラメータを受け取り、パラメータの情報を用いて領域を抽出する（ステップＡ１０）。例えば、ｉ番目のクラスに属する領域を抽出するには、式（１１）で与えられる粗視化されたクラス帰属確率の値をｊ＝１からｊ＝Ｎまで調べ、その確率がある値（閾値）以上の部分空間をｉ番目のクラスに属する部分空間とする。次に、その部分空間に含まれるピクセルを調べ、これらのピクセルをｉ番目の領域であるとして抽出する。クラス帰属確率の閾値としては、例えば0.5を用いれば、所望の結果を得ることができる。
【０１０５】
領域抽出装置８は、本発明請求項４に記載されるように、自動的に閾値を設定することも可能である。そのためには以下のようにする。ｉ番目の領域を抽出するためには、まず推定されたクラス混合比ｗｉに全ピクセル数を掛け、各クラスに属するピクセル数の推定値を求める。この数をｎｉとする。
【０１０６】
次に、式（１１）の粗視化されたクラス帰属確率をｊ＝１からｊ＝Ｎまで調べ、値の高い部分空間から順に、その部分空間に含まれるピクセルを抽出し、抽出したピクセルがｎｉに到達するまで続ける。ｎｉ番目に抽出したピクセルの番号をｌとすると、式（１１）の粗視化されたクラス帰属確率の値が、ｉ番目の領域に属する確率の閾値となる。領域抽出装置８は、領域抽出が完了した後、抽出した領域のデータを出力装置９へ送る。
【０１０７】
上に説明した、データ空間粗視化装置２、粗視化経験分布算出装置３、粗視化条件付き確率分布４、クラス帰属確率算出装置５、パラメータ更新装置６、評価関数算出装置７、領域抽出装置８は、例えば、パーソナルコンピュータやワークステーション、あるいはスーパーコンピュータなどの計算機を用いて構築することができる。出力装置９は、領域抽出装置８から領域データを受け取り、出力する（ステップＡ１１）。出力装置９は、例えばディスプレイ装置やプリンタなどを用いて構成することができる。
【０１０８】
図３は本発明の第２の実施形態の画像領域抽出方法における処理手順を示すフローチャートである。また、図４は本実施形態の画像領域抽出装置の構成を示すブロック図である。本実施形態は、第１の実施形態に領域復元装置１０を追加した点を特徴とする。従って、以下の説明では、第１の実施の形態と重複する部分の説明は省略する。
【０１０９】
領域復元装置１０は、請求項３に記載されるように、粗視化確率分布を用いたパラメータ推定が完了した後、粗視化の解像度が元の解像度に等しいかを調べ（ステップＡ１２）、粗視化されていれば、データを元の解像度に戻す（ステップＡ１３）。粗視化されていなければ、元の解像度での推定が完了していることを意味するから、元の解像度で推定されたパラメータを領域抽出装置８に送る。
【０１１０】
データを元の解像度に戻すには、粗視化の解像度をデータ属性値の最小単位に設定し（例えば１）、第１の実施の形態に述べた方法と全く同一の方法を繰り返せばよい。この場合、第１の実施例で述べた方法よりも推定に時間を要するようになるが、より高い精度でパラメータを推定することができ、結果としてより精度の高い領域抽出が可能となる。また、粗視化確率分布を用いて推定したパラメータは、すでに最適なパラメータの近傍に推定されているため、最初から高い解像度で推定を行うよりも少ない回数のパラメータ更新で最適パラメータが推定でき、はるかに高速に領域抽出を行うことができる。
【０１１１】
次に、シミュレーションデータを用いて本発明の有効性を評価した結果について説明する。
【０１１２】
図８は、細胞の顕微鏡写真を模擬した画像であり、腺腔（画像中もっとも明るい領域）、細胞核（腺腔の周りを取り囲むように存在する輝度の低い小さな円形の組織）、およびそれ以外の組織から成る。書面ではモノクロ画像であるが、実際は腺腔が薄いピンク、細胞核が暗い青紫色、それ以外の領域が赤紫色のカラー画像となっている。
【０１１３】
ここでの目的は、画像から腺腔領域と細胞核領域を分離して抽出することである。抽出した腺腔の形状、細胞核の密度や大きさは、例えば細胞が癌化しているかどうかの判定に役立てることができる。しかしながら図８からわかるように、細胞核に対応する領域は背景と区別が付きにくく、通常の方法では細胞核領域だけを自動抽出することは困難である。
【０１１４】
この画像に対し、本発明に記載する手法を用いて、腺腔領域と細胞核領域を抽出することを試みた。ここでの推定は、粗視化解像度を16×16×16とし、クラスの数を３と仮定して行っている。クラスの数を３としたのは、細胞核に対応する領域（クラス１とする）、腺腔に対応する領域（クラス２とする）、それ以外の領域（クラス３とする）に分離できることを期待しているからである。
【０１１５】
このようにして得られた推定結果を図９（ａ）、（ｂ）、（ｃ）、（ｄ）に示す。図９（ａ）は、クラス１（細胞核領域）のクラス帰属確率を３次元的に表現したものであり、ｘ軸、ｙ軸が画像上のピクセル位置を示し、ｚ軸が帰属確率値を表している。図９（ｂ）はクラス１への帰属確率を２次元的に濃淡で表示したものであり、クラス１に属する確率が低いピクセルほど黒く、確率が高いピクセルほど白く表現されている。
【０１１６】
同様にして図９の（ｃ）、（ｄ）は、クラス２（腺腔領域）へのクラス帰属確率を、それぞれ３次元的、２次元的に表したものであり、帰属確率の高い部分が腺腔領域に属するピクセルとなっていることが分かる。（クラス３に対する結果は、今の場合抽出対象ではないので省略する。）これらの結果は、閾値を陽に設定することなく自動的に得られたものであり、各クラスへの帰属確率の分布を調べるだけで、領域抽出が可能となっている。
【０１１７】
同じ推定を、粗視化無し（粗視化解像度１）で行った結果を図１０の（ａ）、（ｂ）、（ｃ）、（ｄ）に示す。図１０（ａ）、（ｂ）は、クラス１へのクラス帰属確率を、それぞれ３次元的、２次元的に表したもの、また図１０（ｃ）、（ｄ）は、クラス２への帰属確率を、それぞれ３次元的、２次元的に表したものである。これらの図から、粗視化解像度１６での推定結果と粗視化無しでの推定結果がほとんど同じであることが分かる。
【０１１８】
推定されたパラメータ値を比較すると、例えばクラス１に属する領域の平均ＲＧＢ値は、粗視化無しで推定した場合ＲＧＢ＝（１４３.０６８，８６.６８９，１３５.７３７）、解像度１６で推定した場合ＲＧＢ＝（１４１.５２２，８６.１２８，１３５.５１３）となっており、粗視化の影響はほとんどない。
【０１１９】
しかしながら、粗視化のある無しでは推定時間に大きな差が生じる。この違いを調べるために、さまざまな粗視化解像度で推定したときの推定時間を比較したグラフを図１１に示す。このグラフより、粗視化解像度１６で推定した場合、粗視化無しで推定したときに比べて約１／７の時間で推定できていることがわかる。
【０１２０】
【発明の効果】
以上に説明したように本発明では、クラス帰属確率に基づいて領域抽出を行うことにより、閾値を陽に設定することなく、自動的に画像データから所望の領域を抽出することができる。また、粗視化確率分布を用いることにより、従来例よりも更に高速に領域抽出を行うことができる。
【図面の簡単な説明】
【図１】本発明の第１の実施の形態における処理手順を示すフローチャートである。
【図２】本発明の第１の実施の形態における、画像領域抽出装置の構成を示すブロック図である。
【図３】本発明の第２の実施の形態における処理手順を示すフローチャートである。
【図４】本発明の第２の実施の形態における、画像領域抽出装置の構成を示すブロック図である。
【図５】本発明の第１の実施の形態における、パラメータの初期化手順を示すフローチャートである。
【図６】粗視化データ空間と部分空間の一形態を示す図である。
【図７】 本発明請求項７記載の、パラメータの初期化方法の一例を表した図である。
【図８】本発明が適応できる画像の一例である。
【図９】本発明を用いた画像領域抽出を、粗視化解像度１６で行った結果を示しており、（ａ）クラス１（細胞核領域）に属するクラス帰属確率分布の３次元表示、（ｂ）クラス１（細胞核領域）に属するクラス帰属確率分布の２次元表示、（ｃ）クラス２（腺腔領域）に属するクラス帰属確率分布の３次元表示、（ｄ）クラス２（腺腔領域）に属するクラス帰属確率分布の２次元表示をそれぞれ示す。
【図１０】本発明を用いた画像領域抽出を、粗視化無し（粗視化解像度１）で行った結果を示しており、（ａ）クラス１（細胞核領域）に属するクラス帰属確率分布の３次元表示、（ｂ）クラス１（細胞核領域）に属するクラス帰属確率分布の２次元表示、（ｃ）クラス２（腺腔領域）に属するクラス帰属確率分布の３次元表示、（ｄ）クラス２（腺腔領域）に属するクラス帰属確率分布の２次元表示をそれぞれ示す。
【図１１】本発明を用いた画像領域抽出を、粗視化解像度１、２、４、８、１６で行ったときの推定時間を示すグラフである。
【符号の説明】
１ 入力装置
２ データ空間粗視化装置
３ 粗視化経験確率分布算出装置
４ 粗視化条件付き確率分布算出装置
５ クラス帰属確率算出装置
６ パラメータ更新装置
７ 評価関数算出装置
８ 領域抽出装置
９ 出力装置
１０ 領域復元装置
１１ 粗視化データ空間
１２ 部分空間
Ａ１〜Ａ１３ ステップ
Ｂ１〜Ｂ１２ ステップ
Claims (16)
- 領域抽出の対象となる像を含む画像データに対し、画像上の各ピクセルを複数のクラスに分類し、所望の領域を抽出する画像領域抽出方法において、
前記画像の前記各ピクセルが取りうる属性値全体から構成されるデータ空間を、所与の解像度で部分空間に分解し、該各部分空間に属性値を取るピクセルの番号と、該ピクセルの属性値の平均値と、該ピクセルの個数とを保持し、粗視化データ空間を構成する第１の工程と、
前記各部分空間の前記ピクセルの個数を前記画像に含まれる全ピクセルの個数で割り、前記粗視化データ空間上の粗視化経験確率分布を算出する第２の工程と、
前記各クラスの属性を規定するクラスパラメータ、前記クラスの数、および前記クラスの混合比率を初期化する第３の工程と、
前記クラスパラメータに基づいて、前記各部分空間の前記平均値における条件付き確率を算出し、該条件付き確率のデータ空間における分布である粗視化条件付き確率分布を算出する第４の工程と、
前記画像を構成する各ピクセルが各前記クラスに帰属する確率であるクラス帰属確率を、前記クラスの混合比率に前記粗視化条件付き確率分布を掛けることにより算出する第５の工程と、
粗視化対数尤度を評価関数として用い、該評価関数を増加させるように前記クラスパラメータ、および前記クラスの混合比率を更新する第６の工程と、
前記粗視化条件付き確率分布と、前記粗視化経験確率分布と、前記第６の工程で更新した前記クラスパラメータおよび前記クラスの混合比率とを用いて、前記評価関数を算出する第７の工程と、
前記評価関数が所与の終了条件を満たしているかどうかを調べる第８の工程と、
前記評価関数が前記所与の終了条件を満たした後、前記クラスパラメータおよび前記混合比率を保持し、前記クラス帰属確率に基づいて、前記ピクセルが属するクラスを決定し、前記所望の領域を抽出する第９の工程とを有し、
前記評価関数が前記所与の条件を満たすまで、前記第４、前記第５、前記第６、前記第７、および前記第８の工程を繰り返すことを特徴とする画像領域抽出方法。 - 前記第８の工程において、前記評価関数が前記所与の終了条件を満たしたとき、前記粗視化の解像度が元の解像度であるかを調べる第１０の工程と、前記粗視化の解像度が元の解像度でなかったとき前記部分空間の解像度を元の解像度へ復元する第１１の工程とを有し、
前記第３の工程において、前記第９の工程で保持されているクラスパラメータおよび前記クラスの混合比率を初期値として用い、所与の条件が満たされるまで、前記第４、前記第５、前記第６、前記第７、および前記第８の工程を繰り返すことを特徴とする請求項１に記載の画像領域抽出方法。 - 前記第９の工程において、推定された前記混合比率に前記画像を構成する全ピクセルの個数を掛けて前記各クラスに属するピクセルの個数を算出し、前記クラス帰属確率が高い順から前記ピクセルを選択し、前記各クラスに属するピクセルを決定することを特徴とする請求項１及至２のいずれかに記載の画像領域抽出方法。
- 前記第６の工程において、前記評価関数としてＡＩＣを用い、前記第７の工程において評価関数を減少させるようパラメータを変更することを特徴とする請求項１及至３のいずれかに記載の画像領域抽出方法。
- 前記第６の工程において、前記評価関数としてＭＤＬを用い、前記第７の工程において該評価関数を減少させるようパラメータを変更することを特徴とする請求項１及至３のいずれかに記載の画像領域抽出方法。
- 前記第６の工程において、前記評価関数としてストラクチュラル・リスク（ Structural Risk ）を用い、前記第７の工程において該評価関数を減少させるようパラメータを変更することを特徴とする請求項１及至３のいずれかに記載の画像領域抽出方法。
- 前記第３の工程は、
前記各部分空間同士が近傍にあるかどうかを定義する近傍半径と、前記クラスの数を設定する第１のステップと、
前記各部分空間ごとに、各部分空間の代表値を設定する第２のステップと、
分類対象部分空間の集合を設定する第３のステップと、
前記分類対象部分空間の集合の中から、前記粗視化経験確率の最も高い部分空間を選択する第４のステップと、
前記粗視化経験確率の最も高い部分空間の代表値との距離が、前記近傍半径以内である代表値を持つ部分空間全てを、近傍集合として選択する第５のステップと、
すでに分類が完了したクラスに含まれる部分空間の代表値と、前記近傍集合に含まれる部分空間の代表値との間の最短距離が、前記近傍半径よりも大きいかどうかを調べる第６のステップと、
もし、すでに分類が完了したクラスに含まれる部分空間の代表値と、前記近傍集合に含まれる部分空間の代表値との間の最短距離が、前記近傍半径よりも大きければ、前記近傍集合を新規クラスとし、前記分類対象部分空間から前記近傍集合を削除して前記第４のステップ以降を繰り返す第７のステップと、
前記最短距離が前記近傍半径以下であれば、前記近傍集合をすでに分類済みのクラスに追加し、前記分類対象部分空間から、前記近傍集合を削除する第８のステップと、
前記分類対象部分空間が空集合であるかどうかを調べる第９のステップと、
もし、前記分類対象部分空間が空集合でなければ、前記第４のステップ以下を繰り返し、前記分類対象部分空間が空集合であれば、分類が完了したクラスの数が所与の数以上あるかどうかを調べる第１０のステップと、
もし分類が完了したクラスの数が所与の数よりも少なければ、前記近傍半径を縮小して前記第３のステップ以降を繰り返す第１１のステップと、
もし、前記分類対象部分空間が空集合であり、かつ分類済みのクラスの数が所与の数よりも多ければ、各クラス内で前記クラスパラメータを算出し、これをクラスパラメータの初期値とし、また各クラス内に含まれる部分空間の数の比率を前記混合比率の初期値とする第１２のステップと、
からなることを特徴とする請求項１及至６のいずれかに記載の画像領域抽出方法。 - 領域抽出の対象となる像を含む画像データに対し、画像上の各ピクセルを複数のクラスに分類し、所望の領域を抽出する画像領域抽出装置において、
画像データを読み込む入力装置と、
前記画像の前記各ピクセルが取りうる属性値全体から構成されるデータ空間を、所与の解像度で部分空間に分解し、該各部分空間に属性値を取るピクセルの番号と、該ピクセルの属性値の平均値と、該ピクセルの個数とを保持し、粗視化データ空間を構成する領域粗視化装置と、
前記各部分空間のピクセルの個数を前記画像に含まれる全ピクセルの個数で割り、前記粗視化データ空間上の粗視化経験確率分布を算出する粗視化経験確率分布算出装置と、
前記各クラスの属性を規定するクラスパラメータ、前記クラスの数、および前記クラスの混合比率を初期化し、前記クラスパラメータに基づいて、前記各部分空間の前記平均値における条件付き確率を算出し、該条件付き確率のデータ空間における分布である粗視化条件付き確率分布を算出する粗視化条件付き確率分布算出装置と、
前記画像を構成する各ピクセルが各前記クラスに帰属する確率であるクラス帰属確率を、前記混合比率に前記粗視化条件付き確率分布を掛けることにより算出するクラス帰属確率算出装置と、
評価関数として粗視化対数尤度を用いるときには、該評価関数を増加させるように、また評価関数としてＡＩＣ、ＭＤＬ、ストラクチュラル・リスクのいずれかを用いる場合は該評価関数を減少させるように前記パラメータを更新するパラメータ更新装置と、
前記粗視化対数尤度、前記ＡＩＣ、前記ＭＤＬ、前記ストラクチュラル・リスクのいずれかを用いて前記評価関数を算出する評価関数算出装置と、
前記評価関数が所与の終了条件を満たしているかどうかを調べ、前記評価関数が前記所与の終了条件を満たした後、前記クラスパラメータを保持し、前記クラス帰属確率に基づいて、前記各ピクセルが属するクラスを決定し、領域を抽出する領域抽出装置と、
抽出した領域を出力する出力装置と、
を有することを特徴とする画像領域抽出装置。 - 前記評価関数が前記所与の終了条件を満たしていることが確認されたのち、粗視化の解像度が元の解像度であるかを調べ、前記データ空間の解像度を元の解像度へ復元する解像度復元装置を有することを特徴とする請求項８に記載の画像領域抽出装置。
- 領域抽出の対象となる像を含む画像データに対し、画像上の各ピクセルを複数のクラスに分類することによって所望の画像領域を抽出する処理をコンピュータに実行させるためのプログラムであって、
前記画像の前記各ピクセルが取りうる属性値全体から構成されるデータ空間を、所与の解像度で部分空間に分解し、該各部分空間に属性値を取るピクセルの番号と、該ピクセルの属性値の平均値と、該ピクセルの個数とを保持し、粗視化データ空間を構成する第１の処理と、
前記各部分空間のピクセルの個数を前記画像に含まれる全ピクセルの個数で割り、前記粗視化データ空間上の粗視化経験確率分布を算出する第２の処理と、
前記各クラスの属性を規定するクラスパラメータ、前記クラスの数、および前記クラスの混合比率を初期化する第３の処理と、
前記クラスパラメータに基づいて、前記各部分空間の前記平均値における条件付き確率を算出し、該条件付き確率のデータ空間における分布である粗視化条件付き確率分布を算出する第４の処理と、
前記画像を構成する各ピクセルが各前記クラスに帰属する確率であるクラス帰属確率を、前記混合比率に前記粗視化条件付き確率分布を掛けることにより算出する第５の処理と、
粗視化対数尤度を評価関数として用い、該評価関数を増加させるように前記クラスパラメータ、および前記クラスの混合比率を更新する第６の処理と、
前記粗視化条件付き確率分布と、前記粗視化経験確率分布と、前記第６の工程で更新した前記クラスパラメータおよび前記クラスの混合比率とを用いて、前記評価関数を算出する第７の処理と、
前記評価関数が所与の終了条件を満たしているかどうかを調べる第８の処理と、
前記評価関数が前記所与の終了条件を満たした後、前記クラスパラメータおよび前記混合比率を保持し、前記クラス帰属確率に基づいて、前記ピクセルが属するクラスを決定し、前記所望の領域を抽出する第９の処理とを含み、
前記評価関数が前記所与の条件を満たすまで、前記第４、前記第５、前記第６、前記第７、および前記第８の処理を繰り返すことを特徴とする画像領域抽出プログラム。 - 前記第８の処理において、前記評価関数が前記所与の終了条件を満たしたとき、前記粗視化の解像度が元の解像度であるかを調べる第１０の処理と、前記粗視化の解像度が元の解像度でなかったとき前記部分空間の解像度を元の解像度へ復元する第１１の処理を含み、
前記第３の処理において前記第９の処理で保持されているクラスパラメータおよび前記クラスの混合比率を初期値として用い、前記評価関数が前記所与の条件を満たすまで、前記第４、前記第５、前記第６、前記第７、および前記第８の処理を繰り返すことを特徴とする請求項１０に記載の画像領域抽出プログラム。 - 前記第９の処理は、推定された前記混合比率に前記画像を構成する全ピクセルの個数を掛けて前記各クラスに属するピクセルの個数を算出し、前記クラス帰属確率が高い順から前記ピクセルを選択し、前記各クラスに属するピクセルを決定する処理を含むことを特徴とする請求項１０及至１１のいずれかに記載の画像領域抽出プログラム。
- 前記第６の処理は、前記評価関数として前記粗視化対数尤度に代えてＡＩＣを算出する処理であり、前記第７の処理は、該評価関数を減少させるようパラメータを変更する処理であることを特徴とする請求項１０及至１２のいずれかに記載の画像領域抽出プログラム。
- 前記第６の処理は、前記評価関数として前記粗視化対数尤度に代えてＭＤＬを算出する処理であり、前記第７の処理は、該評価関数を減少させるようパラメータを変更する処理であることを特徴とする請求項１０及至１２のいずれかに記載の画像領域抽出プログラム。
- 前記第６の処理は、前記評価関数として前記粗視化対数尤度に代えてストラクチュラル・リスク（ Structural Risk ）を算出する処理であり、前記第７の処理は、該評価関数を減少させるようパラメータを変更する処理であることを特徴とする請求項１０及至１２のいずれかに記載の画像領域抽出プログラム。
- 前記第３の処理は、
前記各部分空間同士が近傍にあるかどうかを定義する近傍半径と、前記クラスの数を設定する第１のステップと、
前記各部分空間ごとに、各部分空間の代表値を設定する第２のステップと、
分類対象部分空間の集合を設定する第３のステップと、
前記分類対象部分空間の集合の中から、前記粗視化経験確率の最も高い部分空間を選択する第４のステップと、
前記粗視化経験確率の最も高い部分空間の代表値との距離が、前記近傍半径以内である代表値を持つ部分空間全てを、近傍集合として選択する第５のステップと、
すでに分類が完了したクラスに含まれる部分空間の代表値と、前記近傍集合に含まれる部分空間の代表値との間の最短距離が、前記近傍半径よりも大きいかどうかを調べる第６のステップと、
もし、すでに分類が完了したクラスに含まれる部分空間の代表値と、前記近傍集合に含まれる部分空間の代表値との間の最短距離が、前記近傍半径よりも大きければ、前記近傍集合を新規クラスとし、前記分類対象部分空間から前記近傍集合を削除して前記第４のステップ以降を繰り返す第７のステップと、
前記最短距離が前記近傍半径以下であれば、前記近傍集合をすでに分類済みのクラスに追加し、前記分類対象部分空間から、前記近傍集合を削除する第８のステップと、
前記分類対象部分空間が空集合であるかどうかを調べる第９のステップと、
もし、前記分類対象部分空間が空集合でなければ、前記第４のステップ以下を繰り返し、前記分類対象部分空間が空集合であれば、分類が完了したクラスの数が所与の数以上あるかどうかを調べる第１０のステップと、
もし分類が完了したクラスの数が所与の数よりも少なければ、前記近傍半径を縮小して前記第３のステップ以降を繰り返す第１１のステップと、
もし、前記分類対象部分空間が空集合であり、かつ分類済みのクラスの数が所与の数よりも多ければ、各クラス内で前記クラスパラメータを算出し、これをクラスパラメータの初期値とし、また各クラス内に含まれる部分空間の数の比率を前記クラスの混合比率の初期値とする第１２のステップと、
を含むことを特徴とする請求項１０及至１５のいずれかに記載の画像領域抽出プログラム。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002108608A JP3649328B2 (ja) | 2002-04-10 | 2002-04-10 | 画像領域抽出方法および装置 |
| US10/409,140 US7136540B2 (en) | 2002-04-10 | 2003-04-09 | Picture region extraction method and device |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002108608A JP3649328B2 (ja) | 2002-04-10 | 2002-04-10 | 画像領域抽出方法および装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2003303344A JP2003303344A (ja) | 2003-10-24 |
| JP3649328B2 true JP3649328B2 (ja) | 2005-05-18 |
Family
ID=28786527
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2002108608A Expired - Fee Related JP3649328B2 (ja) | 2002-04-10 | 2002-04-10 | 画像領域抽出方法および装置 |
Country Status (2)
| Country | Link |
|---|---|
| US (1) | US7136540B2 (ja) |
| JP (1) | JP3649328B2 (ja) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2011121900A1 (ja) * | 2010-03-30 | 2011-10-06 | 日本電気株式会社 | 画像処理装置、画像読取装置、画像処理方法及び画像処理プログラム |
Families Citing this family (21)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7327862B2 (en) | 2001-04-30 | 2008-02-05 | Chase Medical, L.P. | System and method for facilitating cardiac intervention |
| US7526112B2 (en) | 2001-04-30 | 2009-04-28 | Chase Medical, L.P. | System and method for facilitating cardiac intervention |
| US7693563B2 (en) | 2003-01-30 | 2010-04-06 | Chase Medical, LLP | Method for image processing and contour assessment of the heart |
| JP4609322B2 (ja) * | 2003-07-29 | 2011-01-12 | 日本電気株式会社 | 染色体状態の評価方法および評価システム |
| JP2005141601A (ja) * | 2003-11-10 | 2005-06-02 | Nec Corp | モデル選択計算装置，動的モデル選択装置，動的モデル選択方法およびプログラム |
| US7333643B2 (en) * | 2004-01-30 | 2008-02-19 | Chase Medical, L.P. | System and method for facilitating cardiac intervention |
| JP4483334B2 (ja) * | 2004-02-18 | 2010-06-16 | 富士ゼロックス株式会社 | 画像処理装置 |
| US7848566B2 (en) | 2004-10-22 | 2010-12-07 | Carnegie Mellon University | Object recognizer and detector for two-dimensional images using bayesian network based classifier |
| JP5088329B2 (ja) | 2007-02-13 | 2012-12-05 | 日本電気株式会社 | 細胞特徴量算出装置および細胞特徴量算出方法 |
| JP4493679B2 (ja) * | 2007-03-29 | 2010-06-30 | 富士フイルム株式会社 | 対象領域抽出方法および装置ならびにプログラム |
| JP5347272B2 (ja) | 2008-01-18 | 2013-11-20 | 日本電気株式会社 | スポット定量装置、スポット定量方法及びプログラム |
| JP5365011B2 (ja) * | 2008-01-29 | 2013-12-11 | 日本電気株式会社 | 病理診断支援装置、病理診断支援方法、およびプログラム |
| EP2309493B1 (en) * | 2009-09-21 | 2013-08-14 | Google, Inc. | Coding and decoding of source signals using constrained relative entropy quantization |
| JP5906071B2 (ja) * | 2011-12-01 | 2016-04-20 | キヤノン株式会社 | 情報処理方法、情報処理装置、および記憶媒体 |
| CN102609721B (zh) * | 2012-02-01 | 2014-06-04 | 北京师范大学 | 遥感影像的聚类方法 |
| KR101370496B1 (ko) * | 2013-09-26 | 2014-03-06 | 한국건설기술연구원 | 복합매질로 이루어진 시편에 대한 X-ray CT 영상의 최소 단위에 존재하는 각 순수매질의 부피비 측정방법 |
| CN113950278A (zh) * | 2019-07-10 | 2022-01-18 | 莎益博网络系统株式会社 | 图像分析装置及图像分析方法 |
| WO2021178287A1 (en) * | 2020-03-02 | 2021-09-10 | Texas Instruments Incorporated | Projector with phase hologram modulator |
| CN111401470B (zh) * | 2020-03-31 | 2023-03-10 | 西安电子科技大学 | 基于特征空间分布的Fisher特征选择方法 |
| US20210366139A1 (en) * | 2020-05-21 | 2021-11-25 | Samsung Electronics Co., Ltd. | Method and apparatus for generating depth image |
| CN112092675B (zh) * | 2020-08-31 | 2022-03-25 | 长城汽车股份有限公司 | 一种电池热失控预警方法、系统及服务器 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6415048B1 (en) * | 1993-10-12 | 2002-07-02 | Schneider Medical Technologies, Inc. | Compositional analysis system |
| JP3534009B2 (ja) | 1999-09-24 | 2004-06-07 | 日本電気株式会社 | 輪郭抽出方法及び装置 |
| US6424960B1 (en) * | 1999-10-14 | 2002-07-23 | The Salk Institute For Biological Studies | Unsupervised adaptation and classification of multiple classes and sources in blind signal separation |
| US20040125877A1 (en) * | 2000-07-17 | 2004-07-01 | Shin-Fu Chang | Method and system for indexing and content-based adaptive streaming of digital video content |
-
2002
- 2002-04-10 JP JP2002108608A patent/JP3649328B2/ja not_active Expired - Fee Related
-
2003
- 2003-04-09 US US10/409,140 patent/US7136540B2/en not_active Expired - Fee Related
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2011121900A1 (ja) * | 2010-03-30 | 2011-10-06 | 日本電気株式会社 | 画像処理装置、画像読取装置、画像処理方法及び画像処理プログラム |
| JP2011210156A (ja) * | 2010-03-30 | 2011-10-20 | Nec Corp | 画像処理装置、画像読取装置、画像処理方法及び画像処理プログラム |
| US9438768B2 (en) | 2010-03-30 | 2016-09-06 | Nec Corporation | Image processing apparatus, image reading apparatus, image processing method and information storage medium |
Also Published As
| Publication number | Publication date |
|---|---|
| US7136540B2 (en) | 2006-11-14 |
| JP2003303344A (ja) | 2003-10-24 |
| US20030194132A1 (en) | 2003-10-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP3649328B2 (ja) | 画像領域抽出方法および装置 | |
| JP7134303B2 (ja) | 顕微鏡スライド画像のための焦点重み付き機械学習分類器誤り予測 | |
| US8687879B2 (en) | Method and apparatus for generating special-purpose image analysis algorithms | |
| DE102008060789A1 (de) | System und Verfahren zur nicht überwachten Detektion und Gleason-Abstufung für ein Prostatakrebspräparat (Whole-Mount) unter Verwendung von NIR Fluoreszenz | |
| JP2006153742A (ja) | 病理診断支援装置、病理診断支援プログラム、病理診断支援方法、及び病理診断支援システム | |
| KR20230125169A (ko) | 조직의 이미지 처리 방법 및 조직의 이미지 처리 시스템 | |
| CN113570619B (zh) | 基于人工智能的计算机辅助胰腺病理图像诊断系统 | |
| Sasmal et al. | A survey on the utilization of Superpixel image for clustering based image segmentation | |
| Tenbrinck et al. | Image segmentation with arbitrary noise models by solving minimal surface problems | |
| WO2022258718A1 (en) | Processing multimodal images of tissue for medical evaluation | |
| CN116805319A (zh) | 提供训练机器学习分割算法的训练数据集的方法和系统 | |
| D’APICE et al. | Variational model with nonstandard growth conditions for restoration of satellite optical images using synthetic aperture radar | |
| US20230096719A1 (en) | Scalable and high precision context-guided segmentation of histological structures including ducts/glands and lumen, cluster of ducts/glands, and individual nuclei in whole slide images of tissue samples from spatial multi-parameter cellular and sub-cellular imaging platforms | |
| JP4609322B2 (ja) | 染色体状態の評価方法および評価システム | |
| CN104766068A (zh) | 一种多规则融合的随机游走舌像提取方法 | |
| CN110739051B (zh) | 利用鼻息肉病理图片建立嗜酸性粒细胞占比模型的方法 | |
| Piening et al. | Learning from small data sets: Patch‐based regularizers in inverse problems for image reconstruction | |
| De et al. | Brain tumor classification from radiology and histopathology using deep features and graph convolutional network | |
| Paroli et al. | Reversible jump Markov Chain Monte Carlo methods and segmentation algorithms in hidden Markov models | |
| Karabağ et al. | Segmentation And Modelling of Helanuclear Envelope | |
| Desquesnes et al. | PdEs-based morphology on graphs for cytological slides segmentation and clustering | |
| 우동준 | Adaptive Image Processing Based on Superpixel | |
| CN114581417A (zh) | 一种医学实验处理装置、系统、方法、设备以及存储介质 | |
| CN117975122A (zh) | 一种基于动态稀疏化对比学习的医学图像分类方法 | |
| CN115497604A (zh) | 胶囊内镜的图像处理方法、装置、设备及存储介质 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20041112 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20041227 |
|
| 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: 20050127 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20050209 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20080225 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090225 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100225 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100225 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110225 Year of fee payment: 6 |
|
| LAPS | Cancellation because of no payment of annual fees |