JP3901540B2

JP3901540B2 - クラスタリングプログラム

Info

Publication number: JP3901540B2
Application number: JP2002039485A
Authority: JP
Inventors: 隆洋中井; 眞蓼沼
Original assignee: ATR Advanced Telecommunications Research Institute International
Current assignee: ATR Advanced Telecommunications Research Institute International
Priority date: 2002-02-18
Filing date: 2002-02-18
Publication date: 2007-04-04
Anticipated expiration: 2022-02-18
Also published as: JP2003242508A

Description

【０００１】
【産業上の利用分野】
この発明はクラスタリング方法に関し、特にたとえば複数の標本点をクラスタリングする、クラスタリング方法に関する。
【０００２】
【従来の技術】
従来のこの種のクラスタリング方法としては、単純クラスタリング方法がある。これは、「画像工学（増補）−画像のエレクトロニクス−（テレビジョン学会教科書シリーズ１）」（コロナ社、２０００年４月２０日発行）の４．６．１節に記述されている。
【０００３】
具体的に説明すると、ｎ個のベクトル（Ｘ₁,Ｘ₂,…,Ｘ_n）をクラスタリングする場合には、まず、任意のベクトル、たとえば、Ｘ_iをとり、これを第１クラスタＣ₁の中心Ｙ₁（Ｙ₁＝Ｘ_i）とする。
【０００４】
次に、ベクトルＸ_jをとり、Ｙ₁とＸ_jとの距離Ｄ_1,jを求める。ここで、Ｄ_i,j＞Ｔであれば、Ｘ_jを第２のクラスタＣ₂の中心Ｙ₂（Ｙ₂＝Ｘ_j）とする。しかし、Ｄ_i,j≦Ｔであれば、Ｘ_j∈Ｃ₁とする。
【０００５】
続いて、Ｘ_kをとり、Ｙ₁およびＹ₂とＸ_kとの距離Ｄ_1,k，Ｄ_2,kを求める。ここで、Ｄ_1,k＞Ｔであり、かつＤ_2,k＞Ｔであれば、Ｘ_kを第３クラスタＣ₃の中心Ｙ₃（Ｙ₃＝Ｘ_k）とする。しかし、Ｄ_1,k≦ＴまたはＤ_2,k≦Ｔであれば、Ｘ_kの中心との距離の短い方のクラスタに所属するものとする。
【０００６】
このような処理をすべてのベクトルについて実行することにより、クラスタリングが完了する。
【０００７】
【発明が解決しようとする課題】
しかし、従来の方法では、決定したクラスタ数が妥当であるかどうか明確でなく、クラスタリングの結果から人間がその妥当性を判断するため、場合によってはあまり好ましくないクラスタ数であることがあった。このため、その後の解析処理等に支障を来たしていた。
【０００８】
それゆえに、この発明の主たる目的は、クラスタリング結果の利便性を向上できる、クラスタリング方法を提供することである。
【０００９】
【課題を解決するための手段】
第１の発明は、複数の標本点をクラスタリングするクラスタリングプログラムであって、コンピュータを、各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離をすべての標本点のそれぞれについて検出する距離検出手段、距離検出手段によって検出したユークリッド距離の第１総和値を求める第１総和値算出手段、第１総和値算出手段によって求めた第１総和値を各クラスタに渡って総和を取った第２総和値を求める第２総和値算出手段、クラスタ数を最大値から最小値に向けて変化させたときに第２総和値算出手段によって求めた第２総和値が急激に変化する個所を挟む前後数箇所の第２総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する最適クラス数決定手段、および最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラムである。
【００１０】
第２の発明は、複数の標本点をクラスタリングするクラスタリングプログラムであって、コンピュータを、各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離をすべての標本点のそれぞれについて検出する距離検出手段、距離検出手段によって検出したユークリッド距離の第１総和値を求める第１総和値算出手段、第１総和値算出手段によって求めた第１総和値を各クラスタに渡って総和を取った第２総和値を求める第２総和値算出手段、クラスタ数を最大値から最小値に向けてｎ（ｎは自然数）個ずつ減少させるとき、前回のクラスタリング結果において最もクラスタ間距離が短い２つのクラスタを１つのクラスタにまとめる操作をｎ回繰り返すことにより、クラスタ数をｎ個減少させたときのクラスタリングを決定して第２総和値を求める第２総和値再計算手段、第２総和値再計算手段によって求めた第２総和値が急激に変化する個所を挟む前後数箇所の第２総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する最適クラスタリング数決定手段、および最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラムである。
【００１１】
第３の発明は、複数の標本点をクラスタリングするクラスタリングプログラムであって、コンピュータを、各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離をすべての標本点のそれぞれについて検出する距離検出手段、距離検出手段によって検出したユークリッド距離の第１総和値を求める第１総和値算出手段、第１総和値算出手段によって求めた第１総和値を各クラスタに渡って総和を取った第２総和値を求める第２総和値算出手段、クラスタ数を最大値から最小値に向けて変化させるときに各クラスタ数に分けるすべての分け方について第２総和値を求める第２総和値再計算手段、第２総和値再計算手段によって求めた複数の第２総和値の中で最も小さいものを第３総和値として算出する第３総和値算出手段、クラスタ数を最大値から最小値に向けて変化させたときに第３総和値算出手段によって算出した第３総和値が急激に変化する前後いずれかの第３総和値に対応するクラスタ数を最適クラスタ数に決定する最適クラス多数決定手段、および最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラムである。
【００１２】
【作用】
第１の発明では、クラスタリングには、重心法が用いられ、まず、クラスタ内の標本点とクラスタの代表点（重心点）との距離を当該クラスタに属する全標本点に渡って総和を取った第１総和値を求める。この第１総和値は、各クラスタについて求められ、それらの総和を取った第２総和値が求められる。たとえば、クラスタ数を最大値から最小値に向かって変化させたとき、第２総和値が急激に変化する前後いずれかの第２総和値に対応するクラスタ数を最適クラスタ数に決定する。そして、決定された最適クラスタ数で標本点がクラスタリングされる。
【００１３】
第２の発明は、第１の発明とほとんど同じであるが、クラスタ数を最大値から最小値に向けてｎ（ｎは自然数）個ずつ減少させるときに、前回のクラスタ間距離が最も短い２つのクラスタを１つにまとめる操作をｎ回繰り返すことにより、クラスタ数をｎ個減少させたときのクラスタリングについて第２総和値を求める。そして、第２総和値が急減に変化する個所を挟む前後数箇所の第２総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する。つまり、前回のクラスタ数におけるクラスタリング結果を利用できるので、演算処理の負担を軽くできる。
【００１４】
第３の発明もまた、第１の発明とほとんど同じであるが、クラスタ数をｎ個減らす場合には、クラスタ数をｎ個減少させたときのすべてのクラスタリングについての第２総和値を求め、その中の最小値を第３総和値とし、クラスタ数を最大値から最小値に向けて変更する場合に第３総和値が急激に変化した前後いずれかの第３総和値に対応するクラスタ数が最適クラスタ数に決定される。つまり、上述の発明よりも正確な（適切な）クラスタ数を決定することができる。
【００１５】
【発明の効果】
この発明によれば、最適なクラスタ数へのクラスタリングが可能であるため、クラスタリング結果の信頼性が高く、その後の処理における利便性を向上させることができる。
【００１６】
この発明の上述の目的,その他の目的,特徴および利点は、図面を参照して行う以下の実施例の詳細な説明から一層明らかとなろう。
【００１７】
【実施例】
この実施例のクラスタリングシステムは、図示は省略するが、パーソナルコンピュータ或いはワークステーションのようなコンピュータとデータベースとを備える。
【００１８】
なお、データベースは、コンピュータの内部に設けられてもよく、直接或いはインターネット等を介してコンピュータと通信可能に外部接続されてもよい。
【００１９】
データベースは、たとえば航空機や衛星から撮影した写真（航空写真、衛星写真）或いは絵画に対応する様々な画像（自然画像）の画像データが記憶（蓄積）される画像データベースである。また、複数の画像データについてのテクスチャの特徴量（後で詳細に説明する、Ｄ，Ｖ，Ｐ）を取得（算出）した結果である結果データも、対応する画像データと関連付けて記憶される。
【００２０】
図１に示すように、この実施例では、画像のサイズは、２０４８画素×２０４８画素である。また、テクスチャの平坦度を判別する単位（画像領域）は、この実施例では、３２画素×３２画素の大きさであり、これを小領域ａ_iとする。したがって、この実施例では、小領域ａ_iの個数は４０９６である。
【００２１】
たとえば、ユーザから標本点（テクスチャ特徴量）の抽出指示が入力されると、コンピュータ、実際には、コンピュータ内部に設けられたＣＰＵ（図示せず）が図２に示すフロー図に従ってテクスチャの特徴量を（抽出）取得する特徴量取得処理を実行して、上述したような結果データ（標本点）を得る。
【００２２】
図２を参照して、ユーザからの指示に応じて、コンピュータのＣＰＵが処理を開始すると、ステップＳ１では、ユーザの指示に従って特徴量を取得すべき所望の画像に対応する画像データを読み出す。続くステップＳ３では、読み出した画像を小領域ａ_iに分割する。上述したように、小領域ａ_iは３２画素×３２画素の大きさであるため、分割される小領域ａ_iの個数は４０９６である。
【００２３】
続いて、ステップＳ５では、ＣＰＵはカウンタのカウント値を初期化（ｉ＝１）し、後で詳細に説明するように、ステップＳ７ではｉ番目の小領域ａ_iについてのテクスチャ特徴量Ｄ（第１特徴量），Ｖ（第２特徴量）およびＰ（第３特徴量）の計算処理を実行して、ステップＳ９でカウント値ｉが４０９６かどうかを判断する。つまり、すべての小領域ａ_iについてのテクスチャ特徴量を計算したかどうかを判断する。
【００２４】
ステップＳ９で“ＮＯ”であれば、つまりカウント値ｉが４０９６でなければ、すべての小領域ａ_iについてのテクスチャ特徴量を終了していないと判断して、ステップＳ１１でカウント値ｉをインクリメント（ｉ＝ｉ＋１）してからステップＳ７に戻る。一方、ステップＳ９で“ＹＥＳ”であれば、つまりカウント値ｉが４０９６であれば、すべての小領域ａ_iについてテクスチャ特徴量の計算処理を終了したと判断して、ステップＳ１３に進む。ステップＳ１３では、４０９６個の小領域ａ_iについてのテクスチャ特徴量（Ｄ，Ｖ，Ｐ）を画像データに関連付けてデータベースに記憶する。
【００２５】
図３は、図２に示したステップＳ７の小領域ａ_iのテクスチャ特徴量の計算処理を示すフロー図であり、特徴量の計算処理が開始されると、ＣＰＵは、ステップＳ２１で小領域ａ_iの画像データをＲ，Ｇ，ＢデータからＹデータに変換する。たとえば、小領域ａ_i上の位置（ｎ，ｍ）における画素データをｑ_i（ｎ，ｍ）とすると、各画素はＲ，Ｇ，Ｂの３種類の画素値を持つので、数１のように表記することができる。
【００２６】
【数１】

【００２７】
ただし、ｎ，ｍは、それぞれ、１≦ｎ≦３２，１≦ｍ≦３２を満たす自然数である。ここで、位置すなわち座標系（ｎ，ｍ）は、画像に対する一般的な座標系（行列座標系）であり、水平位置ｎは右向きを正の方向とし、垂直位置ｍは下向きを正の方向とする。したがって、小領域ａ_iの画像の左上が（１，１）であり、右下が（３２，３２）である。
【００２８】
以後の処理においては、テクスチャは濃淡画像として扱うため、ＣＰＵは、数２に従って各画素のＲ，Ｇ，ＢデータをＹデータに変換する。
【００２９】
【数２】

【００３０】
また、これ以降の処理では、テクスチャの平均的な明るさ（平均輝度）は問題にしないので、次のステップＳ２３では、Ｙデータから小領域ａ_iの平均輝度を除去する。
【００３１】
まず、数３に従って小領域ａ_iの平均輝度を求める。
【００３２】
【数３】

【００３３】
次に、小領域ａ_iの各画素の輝度値Ｙ_i（ｎ，ｍ）から小領域ａ_iの平均輝度を減算した減算値（減算データ）ＹＴ_i（ｎ，ｍ）を数４に従って求める。
【００３４】
【数４】

【００３５】
ただし、１≦ｎ≦３２,１≦ｍ≦３２
次に、ステップＳ２５では、小領域ａ_iの減算データＹＴ_i（ｎ，ｍ）を２次元離散的フーリエ変換（２次元ＤＦＴ）する。つまり、数５に従って演算し、複素数Ｆ_i（ｕ，ｖ）を算出する。
【００３６】
【数５】

【００３７】
ただし、ｊは１単位虚数であり、ｕ，ｖは、それぞれ、−１６≦ｕ≦１６，−１６≦ｖ≦１６を満たす整数である。この複素数Ｆ_i（ｕ，ｖ）は３２画素×３２画素の減算データＹＴ_i（ｎ，ｍ）の２次元周波数（ｕ，ｖ）における周波数成分を表す。
【００３８】
ここで、ｕの単位は[cpw : cycle per picture width]であり、ｖの単位は[cph : cycle per picture hight]である。
【００３９】
一般的に、或る物理量が複素数で表される場合、その物理量の大きさや強さ（強度）は、その複素数の大きさすなわち絶対値で表すことができる。そこで、ＣＰＵは、続くステップＳ２７で２次元周波数平面上に数５によって得られた複素数Ｆ_i（ｕ，ｖ）の絶対値を取る（算出する）。具体的には、数６に従って複素数Ｆ_i（ｕ，ｖ）の絶対値（「Ｗ_i（ｕ，ｖ）」と表記する。）を取り、絶対値Ｗ_i（ｕ，ｖ）をもって２次元周波数成分Ｆ_i（ｕ，ｖ）の強度とする。
【００４０】
【数６】

【００４１】
ただし、上述したように、ｕ，ｖは、それぞれ、−１６≦ｕ≦１６，−１６≦ｖ≦１６を満たす整数である。
【００４２】
また、座標系（ｎ，ｍ）に対応して、２次元周波数座標系（ｕ，ｖ）についても、水平周波数ｕは右向きを正の方向とし、垂直周波数ｖは下向きを正の方向とする。このような２次元周波数座標系（平面上）に、２次元周波数強度分布Ｗ_i（ｕ，ｖ）を等高線表示した一例が図４のように示される。この図４では、グレースケールによって周波数分布の強度を表しており、黒に近いほど、強度が弱く、白に近いほど、強度が強い。
【００４３】
続いて、ステップＳ２９では、この強度分布Ｗ_i（ｕ，ｖ）を２次元周波数平面上の荷重分布とみなし、その重心位置を決定する（求める）。ただし、図４から分かるように、荷重分布（強度分布）は点（原点）対称であるため、重心位置が常に原点となることを避けるように、この実施例では、周波数平面の右半分の領域（水平周波数ｕが０から＋１６cpw、垂直周波数ｖが−１６から＋１６cphの範囲）についての重心Ｇ_iの位置（ｕ_iG，ｖ_iG）が数７に従って求められる。
【００４４】
【数７】

【００４５】
ここで、ＳＷ_iは小領域ａ_iの周波数強度分布の右半分の領域における総荷重を表し、数８で与えられる。また、重心位置の測定例が図４内に点Ｐで示される。
【００４６】
【数８】

【００４７】
なお、この実施例では、重心位置を周波数平面の右半分の領域で求めたが、これに限定する必要はなく、周波数平面の左半分、上半分或いは下半分の領域で求めてもよい。
【００４８】
また、計算の煩わしさを考慮しなければ、上、下、左、右半分の領域以外の場合で、原点を通る直線で周波数平面を２分割した一方の領域（範囲）内で重心位置を求めるようにしてもよい。
【００４９】
次に、この実施例では、テクスチャの方向は問題にしていないので、ステップＳ３１では、ＣＰＵは重心位置によらず，原点から重心位置までの距離Ｄ_iを小領域ａ_iのテクスチャの第１特徴量（テクスチャ特徴量Ｄ）として数９により算出する。
【００５０】
【数９】

【００５１】
続いて、ステップＳ３３で、小領域a_iのテクスチャの第２特徴量（テクスチャ特徴量Ｖ）として、上述したような周波数平面の右半分の領域の周波数範囲で、周波数強度分布の重心周りの分散ｖ_iを数１０により求める。このとき、周波数範囲の強度分布が、この小領域ａ_iに関する全サンプルであるので、不偏分散を求める必要はない。このため、数１０においては、分母がＳＷ_i-１ではなく、ＳＷ_iとされる。
【００５２】
【数１０】

【００５３】
ここで、小領域のテクスチャに関して、右半分の領域における水平周波数ｕを第１変量、垂直周波数ｖを第２変量とすると、この２つの変量ｕ，ｖを次のように統合して新しい２つの変量、すなわち第１主成分z_i1および第２主成分z_i2に変換することができる。
【００５４】
まず、第１変量ｕの分散（水平周波数ｕに沿って分布する荷重による分散）Ｖ_iuu、第２変量ｖの分散（垂直周波数ｖに沿って分布する荷重による分散）Ｖ_ivv、および第１変量ｕ、第２変量ｖの共分散Ｖ_iuvは数１１〜数１３によってそれぞれ求められる。
【００５５】
【数１１】

【００５６】
【数１２】

【００５７】
【数１３】

【００５８】
ここで、これらの分散値を要素とする分散・共分散行列（実対称行列）Ｖ_iは数１４で示される。
【００５９】
【数１４】

【００６０】
この実対称行列Ｖ_iの固有値λ_i1，λ_i2（ただし、λ_i1≧λ_i2≧０である。）および固有値λ_i1，λ_i2に対応する固有ベクトルｌ_i1，ｌ_i2は次のように求められる。具体的には、固有値λ_i1，λ_i2は、Ｖ_iの固有多項式（λに関する２次方程式）、すなわち数１５の根として求められる。また、固有ベクトルｌ_i1，ｌ_i2は、数１６を解くことにより求められる。
【００６１】
【数１５】

【００６２】
【数１６】

【００６３】
次に、新しい変量である第１主成分z_i1および第２主成分z_i2は、ｕ，ｖの１次結合であり、数１７によってそれぞれ定義される。
【００６４】
【数１７】

【００６５】
ここで、固有ベクトルｌ_i1，ｌ_i2は、それぞれ、数１８のように表している。
【００６６】
【数１８】

【００６７】
以上より、ＣＰＵは、ステップＳ３５で、２次元周波数平面の右半分の領域で第１主成分の寄与率Ｐ_iを小領域ａ_iのテクスチャの第３特徴量（テクスチャ特徴量Ｐ）として求める。すなわち、第１主成分の寄与率Ｐ_iは数１９で算出することができる。
【００６８】
【数１９】

【００６９】
しかし、第１固有値λ_i1は、第１主成分z_i1の分散Ｖ_imに等しく、第２固有値λ_i2は、第２主成分z_i2の分散Ｖ_isに等しいことが知られており、また、この実施例では、全変量の個数は２個であることから、第１主成分ｚ_i1の分散Ｖ_imと第２主成分ｚ_i2の分散Ｖ_isの和は、全分散Ｖ_iすなわち重心周りの分散に等しくなる。このため、寄与率Ｐ_iは、第１主成分ｚ_i1の分散Ｖ_imが全分散Ｖ_iに占める割合であるということもでき、数２０が成立する。
【００７０】
【数２０】

【００７１】
なお、この第１主成分ｚ_i1の座標軸は、重心を通るあらゆる軸の中で、その軸に沿った分散が最も大きくなる軸を意味する。第１主成分軸と第２主成分軸との測定例は、図４においてそれぞれ実線と点線とで示される。
【００７２】
以上のように、Ｓ３１，Ｓ３３およびＳ３５の処理を経て、この小領域ａ_iのテクスチャ特徴量すなわち、距離Ｄ_i、分散Ｖ_i、第１主成分の寄与率Ｐ_iが得られる。これらが、上述したように、標本値としてデータベースに記憶される。
【００７３】
たとえば、標本値としての結果データが最適なクラスタ数にクラスタリングされ、そのクラスタリング結果を用いて、テクスチャ平坦度の判別等のようなその後の処理に役立てることができる。
【００７４】
ただし、上述したテクスチャ特徴量は、個別に得た値であり、すなわちそれぞれ単位が異なるため、少なくともクラスタリングの開始時において、３つの特徴量を正規化する必要がある。具体的には、３種類の特徴量Ｄ_i，Ｖ_i，Ｐ_iのそれぞれについて、その４０９６個の測定データに渡る平均値が０、分散が１になるように正規化が行なわれる。以下においては、正規化された特徴量（正規化特徴量）を、それぞれ、ＺＤ_i，ＺＶ_i，ＺＰ_iと表記することにする。
【００７５】
まず、テクスチャ特徴量Ｄ_iは、数２１に従って正規化される。
【００７６】
【数２１】

【００７７】
ここで、Ｄ_M，σ_Dは、それぞれＤ_iの平均および標準偏差を表し、数２２および数２３によって計算される。
【００７８】
【数２２】

【００７９】
【数２３】

【００８０】
また、テクスチャ特徴量Ｖ_iは、数２４に従って正規化される。
【００８１】
【数２４】

【００８２】
ここで、Ｖ_M，σ_Vは、それぞれＶ_iの平均および標準偏差を表し、数２５および数２６によって計算される。
【００８３】
【数２５】

【００８４】
【数２６】

【００８５】
さらに、テクスチャ特徴量Ｐ_iは、数２７に従って正規化される。
【００８６】
【数２７】

【００８７】
ここで、Ｐ_M，σ_Pは、それぞれＰ_iの平均および標準偏差を表し、数２８および数２９によって計算される。
【００８８】
【数２８】

【００８９】
【数２９】

【００９０】
続いて、クラスタリングについて簡単に説明すると、クラスタ数ａがＮ個（全標本個数と同じ数）の場合に、クラスタ間の距離の決め方としてたとえば重心法に基づいてクラスタリングすると、クラスタリングは一意に決まる。
【００９１】
クラスタ数がa（ただし、ａ＝Ｎ，Ｎ−１，…，２，１）のときのｉ番目のクラスタをＣ（ａ，ｉ）と表記する。ここで、ｉ＝１，２，…，ａである。また、クラスタＣ（ａ,ｉ）に含まれる標本点の個数をｂ（ａ，ｉ）とし、クラスタＣ（ａ，ｉ）に含まれるｋ番目の標本点をＳ（ａ，ｉ，ｋ）と表記する。ここで、ｋ＝１，２，…，ｂ（ａ，ｉ）である。さらに、クラスタＣ（ａ，ｉ）の重心をＧ（ａ，ｉ）とし、重心Ｇ（ａ，ｉ）と標本点Ｓ（ａ，ｉ，ｋ）との距離をｄ（ａ，ｉ，ｋ）と表記すると、クラスタＣ（ａ，ｉ）内の各標本点を重心Ｇ（ａ，ｉ）で代表させたことによるクラスタＣ（ａ，ｉ）についてのクラスタリング誤差は、数３０で示される。
【００９２】
【数３０】

【００９３】
なお、以下においては、簡単のため、このＤ（ａ，ｉ）をｉ番目のクラスタＣ（ａ，ｉ）の誤差と呼ぶことにする。
【００９４】
最初のクラスタ数ａ＝Ｎのときには、各クラスタＣ（Ｎ，ｉ）（ただし、ｉ＝１,２,…,Ｎ）は1個の標本点のみを含むので、各クラスタＣ（Ｎ，ｉ）内の標本点と各クラスタの重心Ｇ（Ｎ，ｉ）は一致する。したがって、ｂ（Ｎ，ｉ）＝１，ｋ＝１，ｄ（Ｎ，ｉ，ｋ）＝０となるので、Ｄ（ａ，ｉ）＝０となる。すなわち、ａ＝Ｎのとき、どのクラスタＣ（ａ，ｉ）についても、クラスタＣ（ａ，ｉ）の誤差は０である。
【００９５】
また、ｉ番目のクラスタＣ（ａ，ｉ）の誤差Ｄ（ａ，ｉ）を（クラスタ数がａのときの）全クラスタ（ｉ＝１，２，…，ａ）について加算した誤差Ｄ（ａ，ｉ）の総和ＳＤ（ａ）は、数３１で算出される。
【００９６】
【数３１】

【００９７】
つまり、ＳＤ（ａ）は、Ｎ個の標本点をａ個のクラスタにクラスタリングしたことによるクラスタリングの誤差を表していると考えられる。以下、この実施例では、このＳＤ（ａ）をa個へのクラスタリングの誤差と呼ぶことにする。最初のクラスタ数ａ＝Ｎのときのクラスタリングについては、上述したように各クラスタの誤差Ｄ（ａ,ｉ）は０であるので、ＳＤ（ａ）＝ＳＤ（Ｎ）＝０となる。すなわちクラスタ数ａ＝Ｎのとき、ａ個へのクラスタリングの誤差は０である。
【００９８】
次に、クラスタ数ａ＝Ｎ−１の場合に、重心法に基づいてＮ個の全標本点をクラスタリングすると、Ｎ個の標本点の内、最も近い２点をまとめて１つのクラスタとし、残りのクラスタはクラスタ数ａ＝Ｎの場合のままとして、クラスタ数ａ＝Ｎ−１個へのクラスタリングを決定する。そして、決定されたクラスタリングについて、上記の諸量を求める処理（数３１を用いた誤差の総和の算出処理）を繰り返せば、最も誤差の少ないクラスタ数ａ＝Ｎ−１個へのクラスタリングの誤差ＳＤ（ａ）＝ＳＤ（Ｎ−１）を得ることができる。
【００９９】
続いて、クラスタ数ａ＝Ｎ−２の場合には、重心法に基づいて、現在のクラスタ数ａ＝Ｎ−１個のクラスタリングがクラスタ数ａ＝Ｎ−２個へのクラスタリングに変更される。この場合、Ｎ−１個の重心の内、最も近い２つの重心を１つのクラスタにまとめるように、クラスタ数ａ＝Ｎ−２個へのクラスタリングを決定する。
【０１００】
なお、このようなクラスタリングの分け方の列挙は、一般に流通している統計解析ソフトウェアを使用するか、自分でプログラムを作成すれば、可能であり、さらに自動化することも容易である。以下、クラスタリングする場合については同様のことが言える。また、一般的に流通している統計処理ソフトウェアとしては、ＳＰＳＳｌｎｃ．社製の「ＳＰＳＳ」やＳｔａｔＳｏｆｔ社製の「ＳＴＡＴＩＳＴＩＣＡ」を用いることができる。
【０１０１】
そして、決定されたクラスタリングで、上記の諸量を求める処理（数３１を用いた誤差の総和の算出処理）を繰り返せば、実用上差し支えのない程度に誤差の少ないクラスタ数ａ＝Ｎ−２個へのクラスタリングの誤差ＳＤ（ａ）＝ＳＤ（Ｎ−２）を得ることができる。
【０１０２】
ただし、Ｎ−２個へのクラスタリングを決める場合には、前回のクラスタリングの結果、すなわちＮ−１個へのクラスタリングの結果を用いずに、Ｎ個の標本をＮ−２個のクラスタに分けるすべての分け方について、上記の手順によってＳＤ（ａ）＝ＳＤ（Ｎ−２）を算出し、その中で最も小さいものを真のＳＤ（ａ）とするようにしてもよい。
【０１０３】
この後者の方法の方がＳＤ（ａ）の評価としてはより正確ではあるが、全標本の個数Ｎが非常に大きい場合には、コンピュータの処理能力にもよるが、現実的な方法ではなくなってしまう。このような場合には、前者の方法を用いて、実用上差し支えのない程度に誤差の少ないＳＤ（ａ）を求めるのがよいと考えられる。
【０１０４】
ここで、クラスタに分けるすべての分け方について具体的に説明することにする。たとえば、５個（Ｎ＝５）の標本点を３個のクラスタに分ける（クラスタリングする）すべての分け方は次のように、場合分けすることができる。ただし、５個の標本点（要素）は１〜５の数字で表すことにする。
【０１０５】
（ｉ）要素の数が３，１，１の３つのクラスタに分ける場合には、クラスタの組み合わせは次のようになる。なお、組み合わせ数は１０（₅Ｃ₃）である。
【０１０６】
（１，２，３），（４），（５）
（１，２，４），（３），（５）
（１，２，５），（３），（４）
（１，３，４），（２），（５）
（１，３，５），（２），（４）
（１，４，５），（２），（３）
（２，３，４），（１），（５）
（２，３，５），（１），（４）
（２，４，５），（１），（３）
（３，４，５），（１），（２）
（ii）要素の数が２，２，１の３つのクラスタに分ける場合には、クラスタの組み合わせは次のようになる。
【０１０７】
（１，２），（３，４），（５）
（１，３），（２，４），（５）
（１，４），（２，３），（５）
（１，２），（３，５），（４）
（１，３），（２，５），（４）
（１，５），（２，３），（４）
（１，２），（４，５），（３）
（１，４），（２，５），（３）
（１，５），（２，４），（３）
（１，３），（４，５），（２）
（１，４），（３，５），（２）
（１，５），（３，４），（２）
（２，３），（４，５），（１）
（２，４），（３，５），（１）
（２，５），（３，４），（１）
以上、（ｉ），（ii）で示すように、５個の標本点を３つのクラスタにクラスタリングする場合には、２５通りのクラスタリング方法がある。このように、分け方は多数存在するため、上述したように、前回の結果を利用してクラスタリングすることにより、コンピュータ（ＣＰＵ）の処理量を軽減してもよい。
【０１０８】
さらに、クラスタ数ａ＝Ｎ−３の場合には、重心法に基づいて、現在のクラスタ数ａ＝Ｎ−２個のクラスタリングがクラスタ数ａ＝Ｎ−３個へのクラスタリングに変更される。この場合、Ｎ−２個の重心の内、最も近い２つの重心を１つのクラスタにまとめるように、クラスタ数ａ＝Ｎ−３個へのクラスタリングを決め、上記の諸量を求める処理（数３１を用いた誤差の総和の算出処理）を繰り返せば、実用上差し支えのない程度に誤差の少ないクラスタ数ａ＝Ｎ−３個へのクラスタリングの誤差ＳＤ（ａ）＝ＳＤ（Ｎ−３）を得ることができる。
【０１０９】
ただし、上述した場合と同様に、Ｎ−３個へのクラスタリングを決める場合には、前回のクラスタリングの結果、すなわちＮ−２個へのクラスタリングの結果を用いずに、Ｎ個の標本をＮ−３個のクラスタに分けるすべての分け方について上記の手順にてＳＤ（ａ）＝ＳＤ（Ｎ−３）を算出し、その中で最も小さいものを真のＳＤ（ａ）とするようにしてもよい。
【０１１０】
このようにして、クラスタリングが繰り返され、クラスタ数ａ＝Ｎ−（Ｎ−２）＝２の場合には、重心法に基づいて、現在のクラスタ数ａ＝３個のクラスタリングが２個のクラスタリングに変更される。この場合、３個の重心の内、最も近い２つの重心を１つのクラスタにまとめるように、クラスタ数ａ＝２個へのクラスタリングを決め、上記の諸量を求める処理（数３１を用いた誤差の総和の算出処理）を繰り返せば、実用上差し支えのない程度に誤差の少ないクラスタ数（ａ＝２）へのクラスタリングの誤差ＳＤ（ａ）＝ＳＤ（２）を得ることができる。
【０１１１】
ただし、２個へのクラスタリングを決める場合には、前回のクラスタリングの結果、つまり、３個へのクラスタリングの結果を用いずに、Ｎ個の標本を２個のクラスタに分けるすべての分け方について、上記の手順によってＳＤ（ａ）＝ＳＤ（２）を算出し、その中で最も小さいものを真のＳＤ（ａ）とするようにしてもよい。
【０１１２】
最後に、ａ＝Ｎ−（Ｎ−１）＝１の場合には、重心法に基づいて、現在のクラスタ数ａ＝２個のクラスタリングをクラスタ数ａ＝１個へのクラスタリングに変える。この場合のクラスタリングは一意に決まる。統合された1つのクラスタの重心は２個の重心の中点となる。このクラスタリングについて、上記の諸量を求める処理（数３１を用いた誤差の総和の算出処理）を繰り返せば、最も誤差の少ないａ＝１個へのクラスタリングの誤差ＳＤ（ａ）＝ＳＤ（１）を得ることができる。
【０１１３】
上述のようにして得たａ個へのクラスタリングの誤差ＳＤ（ａ）（ただし、ａ＝Ｎ，Ｎ−１，…，２，１）を、横軸をクラスタ数ａとしてプロットする。この実施例では、Ｎ＝１８であり、前回のクラスタリングの結果を用いる方法で得たＳＤ（ａ）をプロットした一例を図５に示す。
【０１１４】
図５を参照して、最適なクラスタ数（最適クラスタ数）ａ_optは、クラスタ数ａをＮから１に向かって減少させたときに、急激にＳＤ（ａ）が増大（変化）する個所を含む変化前（直前）のクラスタ数ａ、または、それより大きなクラスタ数ａであって、許せる程度の大きさ（許容範囲）のＳＤ（ａ）を与えるクラスタ数ａに決定する。したがって、図５に示す例では、最適クラスタ数ａ_optは「３」である。
【０１１５】
ただし、この実施例では、許容範囲は最大誤差の２０パーセント未満としてあるが、設計者や使用者等によって自由に変更可能である。また、最適クラスタ数ａ_optは、急激にＳＤ（ａ）が変化する変化後（直後）のクラスタ数ａであって、許容範囲を超えないクラスタ数ａに決定するようにしてもよい。したがって、許容範囲の決め方によっては、数個所の最適クラスタ数ａ_optの候補が存在することとなり、その中から１つを決定する場合がある。
【０１１６】
具体的なクラスタリングの処理は、図６および図７のフロー図によって示される。図６を参照して、ユーザによってクラスタリングの指示が与えられると、ＣＰＵは処理を開始し、ステップＳ４１でクラスタリング指示に応じた４０９６個の標本値をデータベースから読み出す。続くステップＳ４３では、読み出した標本値すなわち（テクスチャ）特徴量（Ｄ，Ｖ，Ｐ）を正規化し、正規化した特徴量（正規化特徴量）ＺＤ，ＺＶおよびＺＰを求める。
【０１１７】
そして、ステップＳ４５では、正規化特徴量ＺＤ，ＺＶおよびＺＰを正規化特徴量空間すなわちＺＤ，ＺＶ，ＺＰ（３次元）空間上にプロットし、図８に示すような散布図を描く。
【０１１８】
以下、この実施例において、４０９６個の小領域ａ_iのテクスチャ特徴量Ｄ_i，Ｖ_i，Ｐ_iをまとめて（Ｄ，Ｖ，Ｐ）と表記することとする。同様に、正規化特徴量ＺＤ_i，ＺＶ_i，ＺＰ_iをまとめて（ＺＤ，ＺＶ，ＺＰ）と表記することにする。
【０１１９】
また、上述したステップＳ４５では、４０９６個の正規化特徴量の測定値（ＺＤ，ＺＶ，ＺＰ）が、ＺＤ，ＺＶ，ＺＰ空間上にプロットされる。ただし、この実施例では、簡単に説明するため、或る画像について、１８個の測定値（ＺＤ，ＺＶ，ＺＰ）がプロットされた例を図８に示してある。
【０１２０】
次にステップＳ４７では、クラスタ数ａを初期化（ａ＝Ｎ）する。続いて、ステップＳ４９では、重心法によりＮ個の全標本値をクラスタリングする。そして、ステップＳ５１でカウンタのカウント値ｉを初期化（ｉ＝１）して、ステップＳ５３でクラスタＣ（ａ，ｉ）の重心Ｇ（ａ，ｉ）を算出する。
【０１２１】
ステップＳ５５では、カウンタのカウント値ｋを初期化（ｋ＝１）し、ステップＳ５７では、クラスタＣ（ａ，ｉ）内の重心Ｇ（ａ，ｉ）と各標本点Ｓ（ａ，ｉ，ｋ）との距離ｄ（ａ，ｉ，ｋ）を算出する。
【０１２２】
次に図７に示すように、次のステップＳ５９では、ｋ＝ｂ（ａ，ｉ）かどうかを判断する。つまり、クラスタＣ（ａ，ｉ）における全標本点との距離ｄ（ａ，ｉ，ｋ）を算出したかどうかを判断する。
【０１２３】
ステップＳ５９で“ＮＯ”であれば、つまりｋ＝ｂ（ａ，ｉ）でなければ、クラスタＣ（ａ，ｉ）における全標本点との距離ｄ（ａ，ｉ，ｋ）を算出していないと判断し、ステップＳ６１でカウント値ｋをインクリメント（ｋ＝ｋ＋１）して、図６に示したステップＳ５７に戻る。
【０１２４】
一方、ステップＳ５９で“ＹＥＳ”であれば、つまりｋ＝ｂ（ａ，ｉ）であれば、クラスタＣ（ａ，ｉ）における全標本点との距離ｄ（ａ，ｉ，ｋ）を算出した判断し、ステップＳ６３でｉ番目のクラスタＣ（ａ，ｉ）の誤差Ｄ（ａ，ｉ）を数３０に従って算出する。
【０１２５】
続くステップＳ６５では、カウント値ｉがクラスタ数ａと等しいかどうかを判断する。このステップＳ６５で“ＮＯ”であれば、つまりカウント値ｉがクラスタ数ａと等しくなければ、ステップＳ６７でカウンタをインクリメント（ｉ＝ｉ＋１）して、図６に示したステップＳ５３に戻る。
【０１２６】
一方、ステップＳ６５で“ＹＥＳ”であれば、つまりカウント値ｉがクラスタ数ａと等しければ、ステップＳ６９でａ個へのクラスタリングの誤差ＳＤ（ａ）を算出する。具体的には、数３１を用いて、ａ個のクラスタにクラスタリングしたことによるクラスタリングの誤差を求める。
【０１２７】
続いて、ステップＳ７１では、クラスタ数ａが１かどうかを判断する。このステップＳ７１で“ＮＯ”であれば、つまりクラスタ数ａが１でなければ、ステップＳ７３でクラスタ数ａをデクリメントし（ａ＝ａ−１）し、重心法に基づいて、今回の（デクリメント後の）クラスタ数ａ（ａ個）のクラスタにクラスタリングしてから、図６に示したステップＳ５１に戻る。
【０１２８】
一方、ステップＳ７１で“ＹＥＳ”であれば、つまりクラスタ数ａが１であれば、ステップＳ７５で誤差ＳＤ(ａ)が急激に増大する直前のクラスタ数ａ、または、それより大きなクラスタ数ａで許容できる大きさ（許容範囲）における誤差ＳＤ（ａ）を与えるａを最適クラスタ数ａ_optに決定してから処理を終了する。
【０１２９】
ただし、この実施例では、誤差ＳＤ(ａ)の許容範囲は、最大誤差（約５５）の１／５（約１０）とし、誤差ＳＤ（ａ）が１０を超えない範囲で最適クラスタ数ａ_optが決定される。つまり、許容範囲は、プログラム（ソフト）の設計者や使用者が任意に設定できる値である。
【０１３０】
この実施例によれば、クラスタ数を最適（妥当）な値に決定してクラスタリングできるので、単なる統計として検証する際の信頼性が高く、その後の処理においてクラスタリング結果が悪影響を与えることはない。つまり、クラスタリング結果の利便性を向上させることができる。
【０１３１】
なお、上述の実施例では、クラスタ数ａを１個ずつ減らすようにしたが、全標本個数Ｎが膨大である場合には、クラスタ数ａを複数個（たとえば、１０個）ずつ減らすようにしてもよい。
【０１３２】
また、上述の実施例では、クラスタリングの手法として、ユークリッド距離に基づく重心法を適用した場合についてのみ説明したが、類似度を表す距離としてはユークリッド距離に限定する必要はない。たとえば、ユークリッド距離の２乗、マハラノビス距離、相関係数等を用いるようにしてもよい。
【０１３３】
さらに、クラスタ間の距離の決め方についても、重心法に限定する必要はなく、最短距離法、最長距離法、群平均法、メディアン法またはウォード法等を用いることも可能である。
【０１３４】
さらにまた、上述の実施例では、テクスチャ特徴量（標本値）を正規化するようにしたが、予め正規化されているデータが標本値である場合には正規化する必要はない。また、標本値はテクスチャ特徴量に限定される必要はなく、他の様々なデータを標本値にすることができる。
【０１３５】
また、上述の実施例では、小領域ａ_iに対する前処理としてＲ，Ｇ，ＢデータからＹデータを求めるようにしたが、Ｒ，Ｇ，Ｂデータから色相Ｈ，彩度Ｓ，明度Ｉ（または、Ｖと表記されることもある。）のデータすなわちＨＳＩ空間のデータに変換し、その結果である明度Ｉのデータを用いるようにしてもよい。
【０１３６】
さらに、上述の実施例では、２次元周波数成分の算出にあたり、２次元離散的フーリエ変換（ＤＦＴ）を用いたが、特に２次元高速フーリエ変換（ＦＦＴ）、或いは、２次元離散的コサイン変換（ＤＣＴ）を用いてもよい。
【図面の簡単な説明】
【図１】この発明の一実施例のクラスタリングで取り扱われる画像を示す図解図である。
【図２】テクスチャ特徴量の抽出処理を示すフロー図である。
【図３】テクスチャ特徴量の計算処理を示すフロー図である。
【図４】テクスチャ特徴量の計算処理において、減算データを２次元ＤＦＴして得られた複素数の絶対値を２次元周波数平面上に等高線表示した図解図である。
【図５】クラスタリング処理のクラスタ数に対する誤差の一例を示すグラフである。
【図６】クラスタリング処理の一部を示すフロー図である。
【図７】クラスタリング処理の他の一部を示すフロー図である。
【図８】テクスチャ特徴量の抽出処理で抽出したテクスチャ特徴量に従ってプロットした一例を示すグラフである。

Claims

複数の標本点をクラスタリングするクラスタリングプログラムであって、
コンピュータを、
各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離を前記すべての標本点のそれぞれについて検出する距離検出手段、
前記距離検出手段によって検出したユークリッド距離の第１総和値を求める第１総和値算出手段、
前記第１総和値算出手段によって求めた第１総和値を各クラスタに渡って総和を取った第２総和値を求める第２総和値算出手段、
クラスタ数を最大値から最小値に向けて変化させたときに前記第２総和値算出手段によって求めた第２総和値が急激に変化する個所を挟む前後数箇所の前記第２総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する最適クラス数決定手段、および
前記最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラム。
複数の標本点をクラスタリングするクラスタリングプログラムであって、
コンピュータを、
各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離を前記すべての標本点のそれぞれについて検出する距離検出手段、
前記距離検出手段によって検出したユークリッド距離の第１総和値を求める第１総和値算出手段、
前記第１総和値算出手段によって求めた第１総和値を各クラスタに渡って総和を取った第２総和値を求める第２総和値算出手段、
クラスタ数を最大値から最小値に向けてｎ（ｎは自然数）個ずつ減少させるとき、前回のクラスタリング結果において最もクラスタ間距離が短い２つのクラスタを１つのクラスタにまとめる操作をｎ回繰り返すことにより、クラスタ数をｎ個減少させたときのクラスタリングを決定して前記第２総和値を求める第２総和値再計算手段、
前記第２総和値再計算手段によって求めた第２総和値が急激に変化する個所を挟む前後数箇所の前記第２総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する最適クラスタリング数決定手段、および
前記最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラム。
複数の標本点をクラスタリングするクラスタリングプログラムであって、
コンピュータを、
各クラスタ内のすべての標本点と当該クラスタ内の重心となる代表点とのユークリッド距離を前記すべての標本点のそれぞれについて検出する距離検出手段、
前記距離検出手段によって検出したユークリッド距離の第１総和値を求める第１総和値算出手段、
前記第１総和値算出手段によって求めた第１総和値を各クラスタに渡って総和を取った第２総和値を求める第２総和値算出手段、
クラスタ数を最大値から最小値に向けて変化させるときに各クラスタ数に分けるすべての分け方について前記第２総和値を求める第２総和値再計算手段、
前記第２総和値再計算手段によって求めた複数の前記第２総和値の中で最も小さいものを第３総和値として算出する第３総和値算出手段、
前記クラスタ数を最大値から最小値に向けて変化させたときに前記第３総和値算出手段によって算出した第３総和値が急激に変化する前後いずれかの前記第３総和値に対応するクラスタ数を前記最適クラスタ数に決定する最適クラス多数決定手段、および
前記最適クラスタ数決定手段によって決定した最適クラスタ数へクラスタリングするクラスタリング手段として機能させる、クラスタリングプログラム。