JP3901540B2 - クラスタリングプログラム - Google Patents
クラスタリングプログラム Download PDFInfo
- Publication number
- JP3901540B2 JP3901540B2 JP2002039485A JP2002039485A JP3901540B2 JP 3901540 B2 JP3901540 B2 JP 3901540B2 JP 2002039485 A JP2002039485 A JP 2002039485A JP 2002039485 A JP2002039485 A JP 2002039485A JP 3901540 B2 JP3901540 B2 JP 3901540B2
- Authority
- JP
- Japan
- Prior art keywords
- clustering
- cluster
- value
- clusters
- sum
- 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
Landscapes
- Complex Calculations (AREA)
- Image Analysis (AREA)
Description
【産業上の利用分野】
この発明はクラスタリング方法に関し、特にたとえば複数の標本点をクラスタリングする、クラスタリング方法に関する。
【0002】
【従来の技術】
従来のこの種のクラスタリング方法としては、単純クラスタリング方法がある。これは、「画像工学(増補)−画像のエレクトロニクス−(テレビジョン学会教科書シリーズ1)」(コロナ社、2000年4月20日発行)の4.6.1節に記述されている。
【0003】
具体的に説明すると、n個のベクトル(X1,X2,…,Xn)をクラスタリングする場合には、まず、任意のベクトル、たとえば、Xiをとり、これを第1クラスタC1の中心Y1(Y1=Xi)とする。
【0004】
次に、ベクトルXjをとり、Y1とXjとの距離D1,jを求める。ここで、Di,j>Tであれば、Xjを第2のクラスタC2の中心Y2(Y2=Xj)とする。しかし、Di,j≦Tであれば、Xj∈C1とする。
【0005】
続いて、Xkをとり、Y1およびY2とXkとの距離D1,k,D2,kを求める。ここで、D1,k>Tであり、かつD2,k>Tであれば、Xkを第3クラスタC3の中心Y3(Y3=Xk)とする。しかし、D1,k≦TまたはD2,k≦Tであれば、Xkの中心との距離の短い方のクラスタに所属するものとする。
【0006】
このような処理をすべてのベクトルについて実行することにより、クラスタリングが完了する。
【0007】
【発明が解決しようとする課題】
しかし、従来の方法では、決定したクラスタ数が妥当であるかどうか明確でなく、クラスタリングの結果から人間がその妥当性を判断するため、場合によってはあまり好ましくないクラスタ数であることがあった。このため、その後の解析処理等に支障を来たしていた。
【0008】
それゆえに、この発明の主たる目的は、クラスタリング結果の利便性を向上できる、クラスタリング方法を提供することである。
【0009】
【課題を解決するための手段】
第1の発明は、複数の標本点をクラスタリングするクラスタリングプログラムであって、コンピュータを、各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離をすべての標本点のそれぞれについて検出する距離検出手段、距離検出手段によって検出したユークリッド距離の第1総和値を求める第1総和値算出手段、第1総和値算出手段によって求めた第1総和値を各クラスタに渡って総和を取った第2総和値を求める第2総和値算出手段、クラスタ数を最大値から最小値に向けて変化させたときに第2総和値算出手段によって求めた第2総和値が急激に変化する個所を挟む前後数箇所の第2総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する最適クラス数決定手段、および最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラムである。
【0010】
第2の発明は、複数の標本点をクラスタリングするクラスタリングプログラムであって、コンピュータを、各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離をすべての標本点のそれぞれについて検出する距離検出手段、距離検出手段によって検出したユークリッド距離の第1総和値を求める第1総和値算出手段、第1総和値算出手段によって求めた第1総和値を各クラスタに渡って総和を取った第2総和値を求める第2総和値算出手段、クラスタ数を最大値から最小値に向けてn(nは自然数)個ずつ減少させるとき、前回のクラスタリング結果において最もクラスタ間距離が短い2つのクラスタを1つのクラスタにまとめる操作をn回繰り返すことにより、クラスタ数をn個減少させたときのクラスタリングを決定して第2総和値を求める第2総和値再計算手段、第2総和値再計算手段によって求めた第2総和値が急激に変化する個所を挟む前後数箇所の第2総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する最適クラスタリング数決定手段、および最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラムである。
【0011】
第3の発明は、複数の標本点をクラスタリングするクラスタリングプログラムであって、コンピュータを、各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離をすべての標本点のそれぞれについて検出する距離検出手段、距離検出手段によって検出したユークリッド距離の第1総和値を求める第1総和値算出手段、第1総和値算出手段によって求めた第1総和値を各クラスタに渡って総和を取った第2総和値を求める第2総和値算出手段、クラスタ数を最大値から最小値に向けて変化させるときに各クラスタ数に分けるすべての分け方について第2総和値を求める第2総和値再計算手段、第2総和値再計算手段によって求めた複数の第2総和値の中で最も小さいものを第3総和値として算出する第3総和値算出手段、クラスタ数を最大値から最小値に向けて変化させたときに第3総和値算出手段によって算出した第3総和値が急激に変化する前後いずれかの第3総和値に対応するクラスタ数を最適クラスタ数に決定する最適クラス多数決定手段、および最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラムである。
【0012】
【作用】
第1の発明では、クラスタリングには、重心法が用いられ、まず、クラスタ内の標本点とクラスタの代表点(重心点)との距離を当該クラスタに属する全標本点に渡って総和を取った第1総和値を求める。この第1総和値は、各クラスタについて求められ、それらの総和を取った第2総和値が求められる。たとえば、クラスタ数を最大値から最小値に向かって変化させたとき、第2総和値が急激に変化する前後いずれかの第2総和値に対応するクラスタ数を最適クラスタ数に決定する。そして、決定された最適クラスタ数で標本点がクラスタリングされる。
【0013】
第2の発明は、第1の発明とほとんど同じであるが、クラスタ数を最大値から最小値に向けてn(nは自然数)個ずつ減少させるときに、前回のクラスタ間距離が最も短い2つのクラスタを1つにまとめる操作をn回繰り返すことにより、クラスタ数をn個減少させたときのクラスタリングについて第2総和値を求める。そして、第2総和値が急減に変化する個所を挟む前後数箇所の第2総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する。つまり、前回のクラスタ数におけるクラスタリング結果を利用できるので、演算処理の負担を軽くできる。
【0014】
第3の発明もまた、第1の発明とほとんど同じであるが、クラスタ数をn個減らす場合には、クラスタ数をn個減少させたときのすべてのクラスタリングについての第2総和値を求め、その中の最小値を第3総和値とし、クラスタ数を最大値から最小値に向けて変更する場合に第3総和値が急激に変化した前後いずれかの第3総和値に対応するクラスタ数が最適クラスタ数に決定される。つまり、上述の発明よりも正確な(適切な)クラスタ数を決定することができる。
【0015】
【発明の効果】
この発明によれば、最適なクラスタ数へのクラスタリングが可能であるため、クラスタリング結果の信頼性が高く、その後の処理における利便性を向上させることができる。
【0016】
この発明の上述の目的,その他の目的,特徴および利点は、図面を参照して行う以下の実施例の詳細な説明から一層明らかとなろう。
【0017】
【実施例】
この実施例のクラスタリングシステムは、図示は省略するが、パーソナルコンピュータ或いはワークステーションのようなコンピュータとデータベースとを備える。
【0018】
なお、データベースは、コンピュータの内部に設けられてもよく、直接或いはインターネット等を介してコンピュータと通信可能に外部接続されてもよい。
【0019】
データベースは、たとえば航空機や衛星から撮影した写真(航空写真、衛星写真)或いは絵画に対応する様々な画像(自然画像)の画像データが記憶(蓄積)される画像データベースである。また、複数の画像データについてのテクスチャの特徴量(後で詳細に説明する、D,V,P)を取得(算出)した結果である結果データも、対応する画像データと関連付けて記憶される。
【0020】
図1に示すように、この実施例では、画像のサイズは、2048画素×2048画素である。また、テクスチャの平坦度を判別する単位(画像領域)は、この実施例では、32画素×32画素の大きさであり、これを小領域aiとする。したがって、この実施例では、小領域aiの個数は4096である。
【0021】
たとえば、ユーザから標本点(テクスチャ特徴量)の抽出指示が入力されると、コンピュータ、実際には、コンピュータ内部に設けられたCPU(図示せず)が図2に示すフロー図に従ってテクスチャの特徴量を(抽出)取得する特徴量取得処理を実行して、上述したような結果データ(標本点)を得る。
【0022】
図2を参照して、ユーザからの指示に応じて、コンピュータのCPUが処理を開始すると、ステップS1では、ユーザの指示に従って特徴量を取得すべき所望の画像に対応する画像データを読み出す。続くステップS3では、読み出した画像を小領域aiに分割する。上述したように、小領域aiは32画素×32画素の大きさであるため、分割される小領域aiの個数は4096である。
【0023】
続いて、ステップS5では、CPUはカウンタのカウント値を初期化(i=1)し、後で詳細に説明するように、ステップS7ではi番目の小領域aiについてのテクスチャ特徴量D(第1特徴量),V(第2特徴量)およびP(第3特徴量)の計算処理を実行して、ステップS9でカウント値iが4096かどうかを判断する。つまり、すべての小領域aiについてのテクスチャ特徴量を計算したかどうかを判断する。
【0024】
ステップS9で“NO”であれば、つまりカウント値iが4096でなければ、すべての小領域aiについてのテクスチャ特徴量を終了していないと判断して、ステップS11でカウント値iをインクリメント(i=i+1)してからステップS7に戻る。一方、ステップS9で“YES”であれば、つまりカウント値iが4096であれば、すべての小領域aiについてテクスチャ特徴量の計算処理を終了したと判断して、ステップS13に進む。ステップS13では、4096個の小領域aiについてのテクスチャ特徴量(D,V,P)を画像データに関連付けてデータベースに記憶する。
【0025】
図3は、図2に示したステップS7の小領域aiのテクスチャ特徴量の計算処理を示すフロー図であり、特徴量の計算処理が開始されると、CPUは、ステップS21で小領域aiの画像データをR,G,BデータからYデータに変換する。たとえば、小領域ai上の位置(n,m)における画素データをqi(n,m)とすると、各画素はR,G,Bの3種類の画素値を持つので、数1のように表記することができる。
【0026】
【数1】
【0027】
ただし、n,mは、それぞれ、1≦n≦32,1≦m≦32を満たす自然数である。ここで、位置すなわち座標系(n,m)は、画像に対する一般的な座標系(行列座標系)であり、水平位置nは右向きを正の方向とし、垂直位置mは下向きを正の方向とする。したがって、小領域aiの画像の左上が(1,1)であり、右下が(32,32)である。
【0028】
以後の処理においては、テクスチャは濃淡画像として扱うため、CPUは、数2に従って各画素のR,G,BデータをYデータに変換する。
【0029】
【数2】
【0030】
また、これ以降の処理では、テクスチャの平均的な明るさ(平均輝度)は問題にしないので、次のステップS23では、Yデータから小領域aiの平均輝度を除去する。
【0031】
まず、数3に従って小領域aiの平均輝度を求める。
【0032】
【数3】
【0033】
次に、小領域aiの各画素の輝度値Yi(n,m)から小領域aiの平均輝度を減算した減算値(減算データ)YTi(n,m)を数4に従って求める。
【0034】
【数4】
【0035】
ただし、1≦n≦32,1≦m≦32
次に、ステップS25では、小領域aiの減算データYTi(n,m)を2次元離散的フーリエ変換(2次元DFT)する。つまり、数5に従って演算し、複素数Fi(u,v)を算出する。
【0036】
【数5】
【0037】
ただし、jは1単位虚数であり、u,vは、それぞれ、−16≦u≦16,−16≦v≦16を満たす整数である。この複素数Fi(u,v)は32画素×32画素の減算データYTi(n,m)の2次元周波数(u,v)における周波数成分を表す。
【0038】
ここで、uの単位は[cpw : cycle per picture width]であり、vの単位は[cph : cycle per picture hight]である。
【0039】
一般的に、或る物理量が複素数で表される場合、その物理量の大きさや強さ(強度)は、その複素数の大きさすなわち絶対値で表すことができる。そこで、CPUは、続くステップS27で2次元周波数平面上に数5によって得られた複素数Fi(u,v)の絶対値を取る(算出する)。具体的には、数6に従って複素数Fi(u,v)の絶対値(「Wi(u,v)」と表記する。)を取り、絶対値Wi(u,v)をもって2次元周波数成分Fi(u,v)の強度とする。
【0040】
【数6】
【0041】
ただし、上述したように、u,vは、それぞれ、−16≦u≦16,−16≦v≦16を満たす整数である。
【0042】
また、座標系(n,m)に対応して、2次元周波数座標系(u,v)についても、水平周波数uは右向きを正の方向とし、垂直周波数vは下向きを正の方向とする。このような2次元周波数座標系(平面上)に、2次元周波数強度分布Wi(u,v)を等高線表示した一例が図4のように示される。この図4では、グレースケールによって周波数分布の強度を表しており、黒に近いほど、強度が弱く、白に近いほど、強度が強い。
【0043】
続いて、ステップS29では、この強度分布Wi(u,v)を2次元周波数平面上の荷重分布とみなし、その重心位置を決定する(求める)。ただし、図4から分かるように、荷重分布(強度分布)は点(原点)対称であるため、重心位置が常に原点となることを避けるように、この実施例では、周波数平面の右半分の領域(水平周波数uが0から+16cpw、垂直周波数vが−16から+16cphの範囲)についての重心Giの位置(uiG,viG)が数7に従って求められる。
【0044】
【数7】
【0045】
ここで、SWiは小領域aiの周波数強度分布の右半分の領域における総荷重を表し、数8で与えられる。また、重心位置の測定例が図4内に点Pで示される。
【0046】
【数8】
【0047】
なお、この実施例では、重心位置を周波数平面の右半分の領域で求めたが、これに限定する必要はなく、周波数平面の左半分、上半分或いは下半分の領域で求めてもよい。
【0048】
また、計算の煩わしさを考慮しなければ、上、下、左、右半分の領域以外の場合で、原点を通る直線で周波数平面を2分割した一方の領域(範囲)内で重心位置を求めるようにしてもよい。
【0049】
次に、この実施例では、テクスチャの方向は問題にしていないので、ステップS31では、CPUは重心位置によらず,原点から重心位置までの距離Diを小領域aiのテクスチャの第1特徴量(テクスチャ特徴量D)として数9により算出する。
【0050】
【数9】
【0051】
続いて、ステップS33で、小領域aiのテクスチャの第2特徴量(テクスチャ特徴量V)として、上述したような周波数平面の右半分の領域の周波数範囲で、周波数強度分布の重心周りの分散viを数10により求める。このとき、周波数範囲の強度分布が、この小領域aiに関する全サンプルであるので、不偏分散を求める必要はない。このため、数10においては、分母がSWi-1ではなく、SWiとされる。
【0052】
【数10】
【0053】
ここで、小領域のテクスチャに関して、右半分の領域における水平周波数uを第1変量、垂直周波数vを第2変量とすると、この2つの変量u,vを次のように統合して新しい2つの変量、すなわち第1主成分zi1および第2主成分zi2に変換することができる。
【0054】
まず、第1変量uの分散(水平周波数uに沿って分布する荷重による分散)Viuu、第2変量vの分散(垂直周波数vに沿って分布する荷重による分散)Vivv、および第1変量u、第2変量vの共分散Viuvは数11〜数13によってそれぞれ求められる。
【0055】
【数11】
【0056】
【数12】
【0057】
【数13】
【0058】
ここで、これらの分散値を要素とする分散・共分散行列(実対称行列)Viは数14で示される。
【0059】
【数14】
【0060】
この実対称行列Viの固有値λi1,λi2(ただし、λi1≧λi2≧0である。)および固有値λi1,λi2に対応する固有ベクトルli1,li2は次のように求められる。具体的には、固有値λi1,λi2は、Viの固有多項式(λに関する2次方程式)、すなわち数15の根として求められる。また、固有ベクトルli1,li2は、数16を解くことにより求められる。
【0061】
【数15】
【0062】
【数16】
【0063】
次に、新しい変量である第1主成分zi1および第2主成分zi2は、u,vの1次結合であり、数17によってそれぞれ定義される。
【0064】
【数17】
【0065】
ここで、固有ベクトルli1,li2は、それぞれ、数18のように表している。
【0066】
【数18】
【0067】
以上より、CPUは、ステップS35で、2次元周波数平面の右半分の領域で第1主成分の寄与率Piを小領域aiのテクスチャの第3特徴量(テクスチャ特徴量P)として求める。すなわち、第1主成分の寄与率Piは数19で算出することができる。
【0068】
【数19】
【0069】
しかし、第1固有値λi1は、第1主成分zi1の分散Vimに等しく、第2固有値λi2は、第2主成分zi2の分散Visに等しいことが知られており、また、この実施例では、全変量の個数は2個であることから、第1主成分zi1の分散Vimと第2主成分zi2の分散Visの和は、全分散Viすなわち重心周りの分散に等しくなる。このため、寄与率Piは、第1主成分zi1の分散Vimが全分散Viに占める割合であるということもでき、数20が成立する。
【0070】
【数20】
【0071】
なお、この第1主成分zi1の座標軸は、重心を通るあらゆる軸の中で、その軸に沿った分散が最も大きくなる軸を意味する。第1主成分軸と第2主成分軸との測定例は、図4においてそれぞれ実線と点線とで示される。
【0072】
以上のように、S31,S33およびS35の処理を経て、この小領域aiのテクスチャ特徴量すなわち、距離Di、分散Vi、第1主成分の寄与率Piが得られる。これらが、上述したように、標本値としてデータベースに記憶される。
【0073】
たとえば、標本値としての結果データが最適なクラスタ数にクラスタリングされ、そのクラスタリング結果を用いて、テクスチャ平坦度の判別等のようなその後の処理に役立てることができる。
【0074】
ただし、上述したテクスチャ特徴量は、個別に得た値であり、すなわちそれぞれ単位が異なるため、少なくともクラスタリングの開始時において、3つの特徴量を正規化する必要がある。具体的には、3種類の特徴量Di,Vi,Piのそれぞれについて、その4096個の測定データに渡る平均値が0、分散が1になるように正規化が行なわれる。以下においては、正規化された特徴量(正規化特徴量)を、それぞれ、ZDi,ZVi,ZPiと表記することにする。
【0075】
まず、テクスチャ特徴量Diは、数21に従って正規化される。
【0076】
【数21】
【0077】
ここで、DM,σDは、それぞれDiの平均および標準偏差を表し、数22および数23によって計算される。
【0078】
【数22】
【0079】
【数23】
【0080】
また、テクスチャ特徴量Viは、数24に従って正規化される。
【0081】
【数24】
【0082】
ここで、VM,σVは、それぞれViの平均および標準偏差を表し、数25および数26によって計算される。
【0083】
【数25】
【0084】
【数26】
【0085】
さらに、テクスチャ特徴量Piは、数27に従って正規化される。
【0086】
【数27】
【0087】
ここで、PM,σPは、それぞれPiの平均および標準偏差を表し、数28および数29によって計算される。
【0088】
【数28】
【0089】
【数29】
【0090】
続いて、クラスタリングについて簡単に説明すると、クラスタ数aがN個(全標本個数と同じ数)の場合に、クラスタ間の距離の決め方としてたとえば重心法に基づいてクラスタリングすると、クラスタリングは一意に決まる。
【0091】
クラスタ数がa(ただし、a=N,N−1,…,2,1)のときのi番目のクラスタをC(a,i)と表記する。ここで、i=1,2,…,aである。また、クラスタC(a,i)に含まれる標本点の個数をb(a,i)とし、クラスタC(a,i)に含まれるk番目の標本点をS(a,i,k)と表記する。ここで、k=1,2,…,b(a,i)である。さらに、クラスタC(a,i)の重心をG(a,i)とし、重心G(a,i)と標本点S(a,i,k)との距離をd(a,i,k)と表記すると、クラスタC(a,i)内の各標本点を重心G(a,i)で代表させたことによるクラスタC(a,i)についてのクラスタリング誤差は、数30で示される。
【0092】
【数30】
【0093】
なお、以下においては、簡単のため、このD(a,i)をi番目のクラスタC(a,i)の誤差と呼ぶことにする。
【0094】
最初のクラスタ数a=Nのときには、各クラスタC(N,i)(ただし、i=1,2,…,N)は1個の標本点のみを含むので、各クラスタC(N,i)内の標本点と各クラスタの重心G(N,i)は一致する。したがって、b(N,i)=1,k=1,d(N,i,k)=0となるので、D(a,i)=0となる。すなわち、a=Nのとき、どのクラスタC(a,i)についても、クラスタC(a,i)の誤差は0である。
【0095】
また、i番目のクラスタC(a,i)の誤差D(a,i)を(クラスタ数がaのときの)全クラスタ(i=1,2,…,a)について加算した誤差D(a,i)の総和SD(a)は、数31で算出される。
【0096】
【数31】
【0097】
つまり、SD(a)は、N個の標本点をa個のクラスタにクラスタリングしたことによるクラスタリングの誤差を表していると考えられる。以下、この実施例では、このSD(a)をa個へのクラスタリングの誤差と呼ぶことにする。最初のクラスタ数a=Nのときのクラスタリングについては、上述したように各クラスタの誤差D(a,i)は0であるので、SD(a)=SD(N)=0となる。すなわちクラスタ数a=Nのとき、a個へのクラスタリングの誤差は0である。
【0098】
次に、クラスタ数a=N−1の場合に、重心法に基づいてN個の全標本点をクラスタリングすると、N個の標本点の内、最も近い2点をまとめて1つのクラスタとし、残りのクラスタはクラスタ数a=Nの場合のままとして、クラスタ数a=N−1個へのクラスタリングを決定する。そして、決定されたクラスタリングについて、上記の諸量を求める処理(数31を用いた誤差の総和の算出処理)を繰り返せば、最も誤差の少ないクラスタ数a=N−1個へのクラスタリングの誤差SD(a)=SD(N−1)を得ることができる。
【0099】
続いて、クラスタ数a=N−2の場合には、重心法に基づいて、現在のクラスタ数a=N−1個のクラスタリングがクラスタ数a=N−2個へのクラスタリングに変更される。この場合、N−1個の重心の内、最も近い2つの重心を1つのクラスタにまとめるように、クラスタ数a=N−2個へのクラスタリングを決定する。
【0100】
なお、このようなクラスタリングの分け方の列挙は、一般に流通している統計解析ソフトウェアを使用するか、自分でプログラムを作成すれば、可能であり、さらに自動化することも容易である。以下、クラスタリングする場合については同様のことが言える。また、一般的に流通している統計処理ソフトウェアとしては、SPSS lnc.社製の「SPSS」やStatSoft社製の「STATISTICA」を用いることができる。
【0101】
そして、決定されたクラスタリングで、上記の諸量を求める処理(数31を用いた誤差の総和の算出処理)を繰り返せば、実用上差し支えのない程度に誤差の少ないクラスタ数a=N−2個へのクラスタリングの誤差SD(a)=SD(N−2)を得ることができる。
【0102】
ただし、N−2個へのクラスタリングを決める場合には、前回のクラスタリングの結果、すなわちN−1個へのクラスタリングの結果を用いずに、N個の標本をN−2個のクラスタに分けるすべての分け方について、上記の手順によってSD(a)=SD(N−2)を算出し、その中で最も小さいものを真のSD(a)とするようにしてもよい。
【0103】
この後者の方法の方がSD(a)の評価としてはより正確ではあるが、全標本の個数Nが非常に大きい場合には、コンピュータの処理能力にもよるが、現実的な方法ではなくなってしまう。このような場合には、前者の方法を用いて、実用上差し支えのない程度に誤差の少ないSD(a)を求めるのがよいと考えられる。
【0104】
ここで、クラスタに分けるすべての分け方について具体的に説明することにする。たとえば、5個(N=5)の標本点を3個のクラスタに分ける(クラスタリングする)すべての分け方は次のように、場合分けすることができる。ただし、5個の標本点(要素)は1〜5の数字で表すことにする。
【0105】
(i)要素の数が3,1,1の3つのクラスタに分ける場合には、クラスタの組み合わせは次のようになる。なお、組み合わせ数は10(5C3)である。
【0106】
(1,2,3),(4),(5)
(1,2,4),(3),(5)
(1,2,5),(3),(4)
(1,3,4),(2),(5)
(1,3,5),(2),(4)
(1,4,5),(2),(3)
(2,3,4),(1),(5)
(2,3,5),(1),(4)
(2,4,5),(1),(3)
(3,4,5),(1),(2)
(ii)要素の数が2,2,1の3つのクラスタに分ける場合には、クラスタの組み合わせは次のようになる。
【0107】
(1,2),(3,4),(5)
(1,3),(2,4),(5)
(1,4),(2,3),(5)
(1,2),(3,5),(4)
(1,3),(2,5),(4)
(1,5),(2,3),(4)
(1,2),(4,5),(3)
(1,4),(2,5),(3)
(1,5),(2,4),(3)
(1,3),(4,5),(2)
(1,4),(3,5),(2)
(1,5),(3,4),(2)
(2,3),(4,5),(1)
(2,4),(3,5),(1)
(2,5),(3,4),(1)
以上、(i),(ii)で示すように、5個の標本点を3つのクラスタにクラスタリングする場合には、25通りのクラスタリング方法がある。このように、分け方は多数存在するため、上述したように、前回の結果を利用してクラスタリングすることにより、コンピュータ(CPU)の処理量を軽減してもよい。
【0108】
さらに、クラスタ数a=N−3の場合には、重心法に基づいて、現在のクラスタ数a=N−2個のクラスタリングがクラスタ数a=N−3個へのクラスタリングに変更される。この場合、N−2個の重心の内、最も近い2つの重心を1つのクラスタにまとめるように、クラスタ数a=N−3個へのクラスタリングを決め、上記の諸量を求める処理(数31を用いた誤差の総和の算出処理)を繰り返せば、実用上差し支えのない程度に誤差の少ないクラスタ数a=N−3個へのクラスタリングの誤差SD(a)=SD(N−3)を得ることができる。
【0109】
ただし、上述した場合と同様に、N−3個へのクラスタリングを決める場合には、前回のクラスタリングの結果、すなわちN−2個へのクラスタリングの結果を用いずに、N個の標本をN−3個のクラスタに分けるすべての分け方について上記の手順にてSD(a)=SD(N−3)を算出し、その中で最も小さいものを真のSD(a)とするようにしてもよい。
【0110】
このようにして、クラスタリングが繰り返され、クラスタ数a=N−(N−2)=2の場合には、重心法に基づいて、現在のクラスタ数a=3個のクラスタリングが2個のクラスタリングに変更される。この場合、3個の重心の内、最も近い2つの重心を1つのクラスタにまとめるように、クラスタ数a=2個へのクラスタリングを決め、上記の諸量を求める処理(数31を用いた誤差の総和の算出処理)を繰り返せば、実用上差し支えのない程度に誤差の少ないクラスタ数(a=2)へのクラスタリングの誤差SD(a)=SD(2)を得ることができる。
【0111】
ただし、2個へのクラスタリングを決める場合には、前回のクラスタリングの結果、つまり、3個へのクラスタリングの結果を用いずに、N個の標本を2個のクラスタに分けるすべての分け方について、上記の手順によってSD(a)=SD(2)を算出し、その中で最も小さいものを真のSD(a)とするようにしてもよい。
【0112】
最後に、a=N−(N−1)=1の場合には、重心法に基づいて、現在のクラスタ数a=2個のクラスタリングをクラスタ数a=1個へのクラスタリングに変える。この場合のクラスタリングは一意に決まる。統合された1つのクラスタの重心は2個の重心の中点となる。このクラスタリングについて、上記の諸量を求める処理(数31を用いた誤差の総和の算出処理)を繰り返せば、最も誤差の少ないa=1個へのクラスタリングの誤差SD(a)=SD(1)を得ることができる。
【0113】
上述のようにして得たa個へのクラスタリングの誤差SD(a)(ただし、a=N,N−1,…,2,1)を、横軸をクラスタ数aとしてプロットする。この実施例では、N=18であり、前回のクラスタリングの結果を用いる方法で得たSD(a)をプロットした一例を図5に示す。
【0114】
図5を参照して、最適なクラスタ数(最適クラスタ数)aoptは、クラスタ数aをNから1に向かって減少させたときに、急激にSD(a)が増大(変化)する個所を含む変化前(直前)のクラスタ数a、または、それより大きなクラスタ数aであって、許せる程度の大きさ(許容範囲)のSD(a)を与えるクラスタ数aに決定する。したがって、図5に示す例では、最適クラスタ数aoptは「3」である。
【0115】
ただし、この実施例では、許容範囲は最大誤差の20パーセント未満としてあるが、設計者や使用者等によって自由に変更可能である。また、最適クラスタ数aoptは、急激にSD(a)が変化する変化後(直後)のクラスタ数aであって、許容範囲を超えないクラスタ数aに決定するようにしてもよい。したがって、許容範囲の決め方によっては、数個所の最適クラスタ数aoptの候補が存在することとなり、その中から1つを決定する場合がある。
【0116】
具体的なクラスタリングの処理は、図6および図7のフロー図によって示される。図6を参照して、ユーザによってクラスタリングの指示が与えられると、CPUは処理を開始し、ステップS41でクラスタリング指示に応じた4096個の標本値をデータベースから読み出す。続くステップS43では、読み出した標本値すなわち(テクスチャ)特徴量(D,V,P)を正規化し、正規化した特徴量(正規化特徴量)ZD,ZVおよびZPを求める。
【0117】
そして、ステップS45では、正規化特徴量ZD,ZVおよびZPを正規化特徴量空間すなわちZD,ZV,ZP(3次元)空間上にプロットし、図8に示すような散布図を描く。
【0118】
以下、この実施例において、4096個の小領域aiのテクスチャ特徴量Di,Vi,Piをまとめて(D,V,P)と表記することとする。同様に、正規化特徴量ZDi,ZVi,ZPiをまとめて(ZD,ZV,ZP)と表記することにする。
【0119】
また、上述したステップS45では、4096個の正規化特徴量の測定値(ZD,ZV,ZP)が、ZD,ZV,ZP空間上にプロットされる。ただし、この実施例では、簡単に説明するため、或る画像について、18個の測定値(ZD,ZV,ZP)がプロットされた例を図8に示してある。
【0120】
次にステップS47では、クラスタ数aを初期化(a=N)する。続いて、ステップS49では、重心法によりN個の全標本値をクラスタリングする。そして、ステップS51でカウンタのカウント値iを初期化(i=1)して、ステップS53でクラスタC(a,i)の重心G(a,i)を算出する。
【0121】
ステップS55では、カウンタのカウント値kを初期化(k=1)し、ステップS57では、クラスタC(a,i)内の重心G(a,i)と各標本点S(a,i,k)との距離d(a,i,k)を算出する。
【0122】
次に図7に示すように、次のステップS59では、k=b(a,i)かどうかを判断する。つまり、クラスタC(a,i)における全標本点との距離d(a,i,k)を算出したかどうかを判断する。
【0123】
ステップS59で“NO”であれば、つまりk=b(a,i)でなければ、クラスタC(a,i)における全標本点との距離d(a,i,k)を算出していないと判断し、ステップS61でカウント値kをインクリメント(k=k+1)して、図6に示したステップS57に戻る。
【0124】
一方、ステップS59で“YES”であれば、つまりk=b(a,i)であれば、クラスタC(a,i)における全標本点との距離d(a,i,k)を算出した判断し、ステップS63でi番目のクラスタC(a,i)の誤差D(a,i)を数30に従って算出する。
【0125】
続くステップS65では、カウント値iがクラスタ数aと等しいかどうかを判断する。このステップS65で“NO”であれば、つまりカウント値iがクラスタ数aと等しくなければ、ステップS67でカウンタをインクリメント(i=i+1)して、図6に示したステップS53に戻る。
【0126】
一方、ステップS65で“YES”であれば、つまりカウント値iがクラスタ数aと等しければ、ステップS69でa個へのクラスタリングの誤差SD(a)を算出する。具体的には、数31を用いて、a個のクラスタにクラスタリングしたことによるクラスタリングの誤差を求める。
【0127】
続いて、ステップS71では、クラスタ数aが1かどうかを判断する。このステップS71で“NO”であれば、つまりクラスタ数aが1でなければ、ステップS73でクラスタ数aをデクリメントし(a=a−1)し、重心法に基づいて、今回の(デクリメント後の)クラスタ数a(a個)のクラスタにクラスタリングしてから、図6に示したステップS51に戻る。
【0128】
一方、ステップS71で“YES”であれば、つまりクラスタ数aが1であれば、ステップS75で誤差SD(a)が急激に増大する直前のクラスタ数a、または、それより大きなクラスタ数aで許容できる大きさ(許容範囲)における誤差SD(a)を与えるaを最適クラスタ数aoptに決定してから処理を終了する。
【0129】
ただし、この実施例では、誤差SD(a)の許容範囲は、最大誤差(約55)の1/5(約10)とし、誤差SD(a)が10を超えない範囲で最適クラスタ数aoptが決定される。つまり、許容範囲は、プログラム(ソフト)の設計者や使用者が任意に設定できる値である。
【0130】
この実施例によれば、クラスタ数を最適(妥当)な値に決定してクラスタリングできるので、単なる統計として検証する際の信頼性が高く、その後の処理においてクラスタリング結果が悪影響を与えることはない。つまり、クラスタリング結果の利便性を向上させることができる。
【0131】
なお、上述の実施例では、クラスタ数aを1個ずつ減らすようにしたが、全標本個数Nが膨大である場合には、クラスタ数aを複数個(たとえば、10個)ずつ減らすようにしてもよい。
【0132】
また、上述の実施例では、クラスタリングの手法として、ユークリッド距離に基づく重心法を適用した場合についてのみ説明したが、類似度を表す距離としてはユークリッド距離に限定する必要はない。たとえば、ユークリッド距離の2乗、マハラノビス距離、相関係数等を用いるようにしてもよい。
【0133】
さらに、クラスタ間の距離の決め方についても、重心法に限定する必要はなく、最短距離法、最長距離法、群平均法、メディアン法またはウォード法等を用いることも可能である。
【0134】
さらにまた、上述の実施例では、テクスチャ特徴量(標本値)を正規化するようにしたが、予め正規化されているデータが標本値である場合には正規化する必要はない。また、標本値はテクスチャ特徴量に限定される必要はなく、他の様々なデータを標本値にすることができる。
【0135】
また、上述の実施例では、小領域aiに対する前処理としてR,G,BデータからYデータを求めるようにしたが、R,G,Bデータから色相H,彩度S,明度I(または、Vと表記されることもある。)のデータすなわちHSI空間のデータに変換し、その結果である明度Iのデータを用いるようにしてもよい。
【0136】
さらに、上述の実施例では、2次元周波数成分の算出にあたり、2次元離散的フーリエ変換(DFT)を用いたが、特に2次元高速フーリエ変換(FFT)、或いは、2次元離散的コサイン変換(DCT)を用いてもよい。
【図面の簡単な説明】
【図1】この発明の一実施例のクラスタリングで取り扱われる画像を示す図解図である。
【図2】テクスチャ特徴量の抽出処理を示すフロー図である。
【図3】テクスチャ特徴量の計算処理を示すフロー図である。
【図4】テクスチャ特徴量の計算処理において、減算データを2次元DFTして得られた複素数の絶対値を2次元周波数平面上に等高線表示した図解図である。
【図5】クラスタリング処理のクラスタ数に対する誤差の一例を示すグラフである。
【図6】クラスタリング処理の一部を示すフロー図である。
【図7】クラスタリング処理の他の一部を示すフロー図である。
【図8】テクスチャ特徴量の抽出処理で抽出したテクスチャ特徴量に従ってプロットした一例を示すグラフである。
Claims (3)
- 複数の標本点をクラスタリングするクラスタリングプログラムであって、
コンピュータを、
各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離を前記すべての標本点のそれぞれについて検出する距離検出手段、
前記距離検出手段によって検出したユークリッド距離の第1総和値を求める第1総和値算出手段、
前記第1総和値算出手段によって求めた第1総和値を各クラスタに渡って総和を取った第2総和値を求める第2総和値算出手段、
クラスタ数を最大値から最小値に向けて変化させたときに前記第2総和値算出手段によって求めた第2総和値が急激に変化する個所を挟む前後数箇所の前記第2総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する最適クラス数決定手段、および
前記最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラム。 - 複数の標本点をクラスタリングするクラスタリングプログラムであって、
コンピュータを、
各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離を前記すべての標本点のそれぞれについて検出する距離検出手段、
前記距離検出手段によって検出したユークリッド距離の第1総和値を求める第1総和値算出手段、
前記第1総和値算出手段によって求めた第1総和値を各クラスタに渡って総和を取った第2総和値を求める第2総和値算出手段、
クラスタ数を最大値から最小値に向けてn(nは自然数)個ずつ減少させるとき、前回のクラスタリング結果において最もクラスタ間距離が短い2つのクラスタを1つのクラスタにまとめる操作をn回繰り返すことにより、クラスタ数をn個減少させたときのクラスタリングを決定して前記第2総和値を求める第2総和値再計算手段、
前記第2総和値再計算手段によって求めた第2総和値が急激に変化する個所を挟む前後数箇所の前記第2総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する最適クラスタリング数決定手段、および
前記最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラム。 - 複数の標本点をクラスタリングするクラスタリングプログラムであって、
コンピュータを、
各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離を前記すべての標本点のそれぞれについて検出する距離検出手段、
前記距離検出手段によって検出したユークリッド距離の第1総和値を求める第1総和値算出手段、
前記第1総和値算出手段によって求めた第1総和値を各クラスタに渡って総和を取った第2総和値を求める第2総和値算出手段、
クラスタ数を最大値から最小値に向けて変化させるときに各クラスタ数に分けるすべての分け方について前記第2総和値を求める第2総和値再計算手段、
前記第2総和値再計算手段によって求めた複数の前記第2総和値の中で最も小さいものを第3総和値として算出する第3総和値算出手段、
前記クラスタ数を最大値から最小値に向けて変化させたときに前記第3総和値算出手段によって算出した第3総和値が急激に変化する前後いずれかの前記第3総和値に対応するクラスタ数を前記最適クラスタ数に決定する最適クラス多数決定手段、および
前記最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラム。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2002039485A JP3901540B2 (ja) | 2002-02-18 | 2002-02-18 | クラスタリングプログラム |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2002039485A JP3901540B2 (ja) | 2002-02-18 | 2002-02-18 | クラスタリングプログラム |
Publications (2)
Publication Number | Publication Date |
---|---|
JP2003242508A JP2003242508A (ja) | 2003-08-29 |
JP3901540B2 true JP3901540B2 (ja) | 2007-04-04 |
Family
ID=27780492
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2002039485A Expired - Fee Related JP3901540B2 (ja) | 2002-02-18 | 2002-02-18 | クラスタリングプログラム |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP3901540B2 (ja) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP4931220B2 (ja) | 2007-03-12 | 2012-05-16 | インターナショナル・ビジネス・マシーンズ・コーポレーション | 検出装置、システム、プログラムおよび検出方法 |
US10262222B2 (en) * | 2016-04-13 | 2019-04-16 | Sick Inc. | Method and system for measuring dimensions of a target object |
-
2002
- 2002-02-18 JP JP2002039485A patent/JP3901540B2/ja not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
JP2003242508A (ja) | 2003-08-29 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108304820B (zh) | 一种人脸检测方法、装置及终端设备 | |
KR101833953B1 (ko) | 이미지들을 비교하는 방법 및 시스템 | |
CN116188805B (zh) | 海量图像的图像内容分析方法、装置和图像信息网络 | |
WO2006129791A1 (ja) | 画像処理システム、3次元形状推定システム、物***置姿勢推定システム及び画像生成システム | |
CN113420640B (zh) | 红树林高光谱图像分类方法、装置、电子设备及存储介质 | |
US20100266212A1 (en) | Estimating Vanishing Points in Images | |
CN109376689B (zh) | 人群分析方法及装置 | |
CN112668461B (zh) | 一种具有野生动物识别的智能监管*** | |
JPWO2019106850A1 (ja) | Sar画像解析システム、画像処理装置、画像処理方法および画像処理プログラム | |
Haris et al. | Superresolution for UAV images via adaptive multiple sparse representation and its application to 3-D reconstruction | |
CN110793441B (zh) | 一种高精度物体几何尺寸测量方法及装置 | |
JP3901540B2 (ja) | クラスタリングプログラム | |
CN112017221B (zh) | 基于尺度空间的多模态图像配准方法、装置和设备 | |
CN112131968A (zh) | 基于dcnn的双时相遥感影像变化检测方法 | |
CN117671031A (zh) | 双目相机标定方法、装置、设备及存储介质 | |
JP3894802B2 (ja) | テクスチャ平坦度判別プログラム | |
CN110472085A (zh) | 三维图像搜索方法、***、计算机设备和存储介质 | |
CN112380967B (zh) | 一种联合图像信息的空间人造目标光谱解混方法及*** | |
CN106296688B (zh) | 基于全局估计的影像模糊检测方法及*** | |
JP4434868B2 (ja) | 画像分割処理システム | |
JP3483113B2 (ja) | 時系列画像検索方法、装置、および時系列画像検索プログラムを記録した記録媒体 | |
CN113506266A (zh) | 舌头腻苔的检测方法、装置、设备及存储介质 | |
Araújo et al. | Comparing the use of sum and difference histograms and gray levels occurrence matrix for texture descriptors | |
CN111860153A (zh) | 尺度自适应的高光谱图像分类方法及*** | |
US20110091106A1 (en) | Image Processing Method And System |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20060914 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20060926 |
|
A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20061003 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20061107 |
|
A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20061120 |
|
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: 20061226 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20061226 |
|
R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
LAPS | Cancellation because of no payment of annual fees |