(57)【特許請求の範囲】
【請求項１】
テクスチャ平坦度判別プログラムであって、
コンピュータを、
画像の小領域に対応する画像データを２次元フーリエ変換し、得られた複素数の絶対値の分布によって２次元周波数平面上の周波数強度分布として決定する分布決定手段、
前記分布決定手段によって決定された周波数強度分布を荷重分布とみなすことにより、前記２次元周波数平面内における所定領域の前記荷重分布から重心を求める重心算出手段、
前記２次元周波数平面の原点から前記重心算出手段によって求められた重心までの距離で表される第１特徴量を取得する第１特徴量取得手段、および
少なくとも前記第１特徴量取得手段によって取得された第１特徴量に基づいて前記小領域のテクスチャが平坦な模様であるかどうかを判別する判別手段として機能させる、テクスチャ平坦度判別プログラム。
【請求項２】
前記所定領域は前記２次元周波数平面内の上半分、下半分、左半分或いは右半分に相当する領域である、請求項１記載のテクスチャ平坦度判別プログラム。
【請求項３】
前記コンピュータを、前記所定領域における前記重心周りの分散で表される第２特徴量を取得する第２特徴量取得手段および前記所定領域における水平周波数を表す変数と垂直周波数を表す変数とを統合して得られる第１主成分の寄与率で表される第３特徴量を取得する第３特徴量取得手段としてさらに機能させ、
前記判別手段は、前記第１特徴量、前記第２特徴量および前記第３特徴量に基づいて前記小領域のテクスチャが平坦な模様であるかどうかを判別する、請求項２記載のテクスチャ平坦度判別プログラム。
【請求項４】
前記コンピュータを、前記第１特徴量、前記第２特徴量および前記第３特徴量をそれぞれ平均値が０、分散が１の変数に正規化する正規化手段としてさらに機能させ、
前記判別手段は、前記正規化手段によって正規化された第１，第２，第３特徴量に基づいて前記小領域のテクスチャが平坦な模様であるかどうかを判別する、請求項３記載のテクスチャ平坦度判別プログラム。
【請求項５】
前記コンピュータを、前記正規化手段によって正規化された第１，第２，第３特徴量を３次元空間上にプロットする描画手段および複数の前記小領域に対応する前記正規化された第１，第２，第３特徴量が示す標本点をクラスタリングするクラスタリング手段としてさらに機能させる、請求項４記載のテクスチャ平坦度判別プログラム。
【請求項６】
前記クラスタリング手段は、各クラスタ内の前記標本点と前記クラスタの代表点との距離を前記クラスタに属する全標本点に渡って総和をとった第１総和値を求める第１総和値算出手段、前記第１総和値算出手段によって求められた第１総和値を前記各クラスタに渡って総和を取った第２総和値を求める第２総和値算出手段、およびクラスタ数を最大値から最小値に向けて変化させたときの各クラスタ数における前記第２総和値が急激に変化する個所を挟む前後数個所の前記第２総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する最適クラスタ数決定手段を含む、請求項５記載のテクスチャ平坦度判別プログラム。
【請求項７】
前記クラスタリング手段は、各クラスタ内の前記標本点と前記クラスタの代表点との距離を前記クラスタに属する全標本点に渡って総和を取った第１総和値を求める第１総和値算出手段、前記第１総和値算出手段によって求められた第１総和値を前記各クラスタに渡って総和を取った第２総和値を求める第２総和値算出手段、クラスタ数を最大値から最小値に向けてｎ（ｎは自然数）個ずつ減少させるとき、前回のクラスタリング結果において最もクラスタ間距離が短い２つのクラスタを１つのクラスタにまとめる操作をｎ回繰り返すことにより、クラスタ数をｎ個減少させたときの第２総和値を求める第２総和値再計算手段、および前記第２総和値再計算手段によって求められた第２総和値が急激に変化する個所を挟む前後数個所の前記第２総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する最適クラスタ数決定手段を含む、請求項５記載のテクスチャ平坦度判別プログラム。
【請求項８】
前記クラスタリング手段は、各クラスタ内の前記標本点と前記クラスタの代表点との距離を前記クラスタに属する全標本点に渡って総和を取った第１総和値を求める第１総和値算出手段、前記第１総和値算出手段によって求められた第１総和値を前記各クラスタに渡って総和を取った第２総和値を求める第２総和値算出手段、クラスタ数を最大値から最小値に向けて変化させるときに各クラスタ数に分けるすべての分け方について第２総和値を求める第２総和値再計算手段、前記第２総和値再計算手段によって求められた複数の第２総和値の中で最も小さいものを第３総和値とする第３総和値算出手段、および前記クラスタ数を最大値から最小値に向けて変化させたときに前記第３総和値が急激に変化する個所を含む前後数箇所の前記第３総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する最適クラスタ数決定手段を含む、請求項５記載のテクスチャ平坦度判別プログラム。
【請求項９】
前記小領域の画像データは平均値が０である輝度信号である、請求項１ないし８のいずれかに記載のテクスチャ平坦度判別プログラム。
【請求項１０】
前記小領域の画像データは平均値が０である明度信号である、請求項１ないし８のいずれかに記載のテクスチャ平坦度判別プログラム。
【請求項１１】
前記２次元フーリエ変換は２次元ＦＦＴを含む、請求項１ないし１０のいずれかに記載のテクスチャ平坦度判別プログラム。
【請求項１２】
前記２次元フーリエ変換は２次元ＤＦＴを含む、請求項１ないし１０のいずれかに記載のテクスチャ平坦度判別プログラム。
【請求項１３】
前記２次元フーリエ変換に代えて２次元ＤＣＴする、請求項１ないし１０のいずれかに記載のテクスチャ平坦度判別プログラム。
【発明の詳細な説明】
【０００１】
【産業上の利用分野】
この発明はテクスチャ平坦度判別方法に関し、特にたとえば画像上の或る領域のテクスチャが平坦的な模様を持っているかどうかを判別するテクスチャ平坦度判別方法に関する。
【０００２】
【従来の技術】
従来からテクスチャの粗さに関する特徴抽出については、いくつかの方法が知られており、たとえば、「画像工学（増補）−画像のエレクトロニクス−（テレビジョン学会教科書シリーズ１）」（コロナ社、２０００年４月２０日発行）の４．２．８節「テクスチャ特徴」に、その代表的な方法が記述されている。
【０００３】
たとえば、空間自己相関関数を利用する方法（方法１）が記述されている。具体的には、画像の濃度値（「輝度値」ともいう。）をｆ（ｎ，ｍ）としたとき、数１によって空間自己相関関数Ａ（ξ，η；ｊ，ｋ）が求められる。
【０００４】
【数１】
0003894802.tif 000002
【０００５】
この数１は、点（ｊ，ｋ）を中心とする（２ｗ＋１）×（２ｗ＋１）の窓の範囲内での相関を示す。たとえば、テクスチャ領域がξ´，η´の周期性を有しているならば、空間自己相関関数Ａ（ξ´，η´；ｊ，ｋ）は大きな値となる。したがって、空間自己相関関数Ａ（ξ´，η´；ｊ，ｋ）が大きな値となるξ´，η´の絶対値の大小を調べることにより、テクスチャの粗さを知ることができる。
【０００６】
また、フーリエ変換係数による方法（方法２）が記述されている。具体的には、画像ｆ（ｘ，y）のフーリエ変換Ｆ（ｕ，ｖ）は数２によって求められる。
【０００７】
【数２】
0003894802.tif 000003
【０００８】
また、Ｆ（ｕ，ｖ）を用いて、エネルギー密度スペクトルＭ（ｕ，ｖ）は数３によって求められる。
【０００９】
【数３】
0003894802.tif 000004
【００１０】
さらに、このエネルギー密度スペクトルＭ（ｕ，ｖ）は、数４に示すように、極座標系（ρ，θ）で表示することができる。
【００１１】
【数４】
0003894802.tif 000005
【００１２】
また、数５によって、環内エネルギーと呼ばれる物理量Ｓを算出する。
【００１３】
【数５】
0003894802.tif 000006
【００１４】
この物理量Ｓは、テクスチャが全体に細かい場合には、原点から離れたところで大きくなり、逆に原点近傍で小さくなるので、テクスチャ全体の細かさを見る（知る）際の特徴量として利用することができる。
【００１５】
なお、上記文献には、他にもいくつかのテクスチャの粗さに関する特徴抽出方法が記述されるが、ここでは省略することにする。
【００１６】
【発明が解決しようとする課題】
しかし、従来の方法（方法１および方法２）は、基本的にはテクスチャの粗さや細かさを識別するためのものであり、テクスチャの明確さを識別することはできなかった。ここで、テクスチャの明確さの判別について説明すると、たとえばテクスチャが非常に細かくてもそのテクスチャが非常に薄ければ、そのテクスチャを細かいと識別するのではなく、模様が平坦なテクスチャと識別し、逆にテクスチャが非常に粗くてもそのテクスチャが非常にはっきりしていれば、そのテクスチャを粗いと識別するのではなく、模様が明確なテクスチャと識別する。
【００１７】
ただし、上述の方法を次のように応用すれば、或る程度テクスチャの明確さを識別できることは、一般的に知られている。
【００１８】
たとえば、上述の方法１の場合には、自己相関関数を求める際に、数１に示したように、分子を数６で表される画像のパワーで割る（除算する）ことにより、正規化しているので、模様の濃い、薄い（濃さ）の情報が除去される。このため、テクスチャの粗さを識別することはできるが、テクスチャの明確さを識別することはできない。しかし、数７に示すように、正規化しない場合には、模様の濃さの情報が除去されていないので、或る程度テクスチャの明確さを識別することは可能である。
【００１９】
【数６】
0003894802.tif 000007
【００２０】
【数７】
0003894802.tif 000008
【００２１】
また、同様に上述の方法２も基本的にはテクスチャの粗さや細かさを識別するためのものであるが、テクスチャが非常に薄ければ、どの位置ρmについても物理量Ｓは小さいので、或る程度のテクスチャの明確さを識別することは可能である。
【００２２】
しかし、このような方法１および方法２では、テクスチャが幾何学的で規則正しい模様である場合には、或る程度の精度でテクスチャの明確さを識別することが可能であるが、実在の風景を写した風景写真、地表を写した航空写真或いは画家によって描かれた絵画などのように、様々な周波数成分と様々なノイズ成分とを含んでいる一般的な画像（自然画像）に含まれるテクスチャに対しては、安定して正確な判別をすることはできないという問題があった。
【００２３】
たとえば、航空機や衛星から撮影された地表の一枚の写真から適当な大きさの領域を切出して、その土地の利用状況を調べるために、その領域のテクスチャが平坦な模様を持っているか、明確な模様を持っているか、或いはその中間的な模様を持っているかの判断が必要になることがある。このような判断は、人間にとっては容易であるが、対象となる写真や切出した領域の数が膨大である場合には、現実には機械を用いた自動判別が必要となる。
【００２４】
また、視覚機能を持つロボットにより撮像された画像から、テクスチャの平坦度を判別し、ロボット周辺の環境を把握して、たとえば壁や段差などの障害（物）を検知して、ロボット自身の進路等を決定したい場合もある。
【００２５】
それゆえに、この発明の主たる目的は、模様が細かいか粗いかに拘わらず、画像上の或る領域のテクスチャが平坦的な模様であるかどうかを精度よく判別できる、テクスチャ平坦度判別方法を提供することである。
【００２６】
【課題を解決するための手段】
この発明は、テクスチャ平坦度判別プログラムであって、コンピュータを、画像の小領域に対応する画像データを２次元フーリエ変換し、得られた複素数の絶対値の分布によって２次元周波数平面上の周波数強度分布として決定する分布決定手段、分布決定手段によって決定された周波数強度分布を荷重分布とみなすことにより、２次元周波数平面内における所定領域の荷重分布から重心を求める重心算出手段、２次元周波数平面の原点から重心算出手段によって求められた重心までの距離で表される第１特徴量を取得する第１特徴量取得手段、および少なくとも第１特徴量取得手段によって取得された第１特徴量に基づいて小領域のテクスチャが平坦な模様であるかどうかを判別する判別手段として機能させる、テクスチャ平坦度判別プログラムである。
【００２７】
【作用】
画像は所定の大きさの小領域に分割され、その小領域に対応する画像データが、たとえば、２次元ＦＦＴ，２次元ＤＦＴのような２次元フーリエ変換や２次元ＤＣＴによって複素数に変換される。このような複素数の絶対値は、２次元周波数平面上に分布する。この分布を周波数強度分布とし、さらに周波数強度分布を荷重分布とみなした場合に、２次元周波数平面の上半分、下半分、左半分或いは右半分の領域（所定領域）についての荷重分布から重心が求められる。そして、２次元周波数平面の原点と重心との距離で表される第１特徴量が取得され、この第１特徴量に基づいて小領域のテクスチャが平坦な模様であるかどうかが判別される。具体的には、第１特徴量が小さいほど、明確な模様のテクスチャであると判別され、第１特徴量が大きいほど、不明確な（平坦な）模様のテクスチャであると判別される。
【００２８】
たとえば、２次元周波数平面の所定領域の重心周りの分散で表される第２特徴量および所定領域における水平周波数の変数と垂直周波数の変数との一次結合で表される第１主成分の寄与率で表される第３特徴量をさらに取得すれば、第１〜第３特徴量に基づいてテクスチャの平坦度を判別できるので、より正確に判別することができる。
【００２９】
また、第１〜第３特徴量に基づいて平坦度を判別する場合には、それぞれの値を正規化することにより、スケールが揃えられる。具体的には、第１〜第３特徴量は、それぞれの平均値が０で、分散が１の変数に正規化される。
【００３０】
各小領域についての正規化された第１〜第３特徴量を、３次元空間にプロットすると、小領域の特徴量の分布を容易に知ることができる。また、プロットされた標本点（小領域についての正規化された特徴量）は、所定のクラスタ数または最適なクラスタ数にクラスタリングされ、統計の分析等に利用される。
【００３１】
このような第１〜第３特徴量の算出では、画像データの色成分は必要ないため、たとえば、Ｒ，Ｇ，ＢデータがＹデータに変換される。つまり、画像の輝度のみに基づいて第１〜第３特徴量が求められる。ただし、輝度に代えて、明度（Ｉデータ）を用いるようにしてもよい。この場合には、Ｒ，Ｇ，ＢデータがＨＳＩ空間のデータ（Ｈ，Ｓ，Ｉデータ）に変換される。
【００３２】
クラスタリングにあたって、クラスタ間の距離の決め方としてはたとえば重心法が用いられ、まず、クラスタ内の標本点とクラスタの代表点（重心点）との距離を当該クラスタに属する全標本点に渡って総和を取った第１総和値を求める。この第１総和値は、各クラスタについて求められ、それらの総和を取った第２総和値が求められる。たとえば、クラスタ数を最大値から最小値に向かって変化させたとき、第２総和値が急激に変化する個所を含む前後数箇所の第２総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定する。そして、決定された最適クラスタ数で標本点がクラスタリングされる。したがって、クラスタリングの結果から小領域画像についてのテクスチャの平坦度を判別することができる。
【００３３】
また、クラスタ数を最大値から最小値に向けてｎ（ｎは自然数）個ずつ減少させるときに、前回のクラスタ結果において最もクラスタ間距離が短い２つのクラスタを１つのクラスタにまとめる操作をｎ回繰り返すことにより、クラスタ数をｎ個減少させたときの第２総和値を求め、その第２総和値が急激に変化する個所を挟む前後数箇所の第２総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定するようにしてもよい。つまり、前回のクラスタ数の結果を利用できるため、演算処理の負担を軽くすることができる。
【００３４】
或いは、各クラスタ数に分けるすべての分け方についての第２総和値を求め、その中の最小値を第３総和値とし、クラスタ数を最大値から最小値に向けて変更する場合に第３総和値が急激に変化した個所を含む前後数箇所の第３総和値のいずれかに対応するクラスタ数を最適クラスタ数に決定するようにすれば、より正確なクラスタ数を決定することができる。
【００３５】
【発明の効果】
この発明によれば、小領域内の全画素に渡って抽出された特徴量を求めるので、テクスチャの平坦度を正確に判別することができる。
【００３６】
この発明の上述の目的,その他の目的,特徴および利点は、図面を参照して行う以下の実施例の詳細な説明から一層明らかとなろう。
【００３７】
【実施例】
この実施例のテクスチャ平坦度判別システムは、図示は省略するが、パーソナルコンピュータ或いはワークステーションのようなコンピュータとデータベースとを備える。
【００３８】
なお、データベースは、コンピュータの内部に設けられてもよく、直接或いはインターネット等を介してコンピュータと通信可能に外部接続されてもよい。
【００３９】
データベースは、たとえば航空機や衛星から撮影した写真（航空写真、衛星写真）或いは絵画に対応する様々な画像（自然画像）の画像データが記憶（蓄積）される画像データベースである。たとえば、図１に示すように、この実施例では、画像のサイズは、２０４８画素×２０４８画素である。また、テクスチャの平坦度を判別する単位（画像領域）は、この実施例では、３２画素×３２画素の大きさであり、これを小領域ａiとする。したがって、この実施例では、小領域ａiの個数は４０９６である。
【００４０】
たとえば、ユーザからの指示が入力されると、コンピュータ、実際には、コンピュータ内部に設けられたＣＰＵ（図示せず）が図２に示すフロー図に従ってテクスチャの判別処理を実行する。
【００４１】
図２を参照して、ユーザからの指示が入力されると、コンピュータのＣＰＵは処理を開始し、ステップＳ１で、ユーザの指示に従ってテクスチャを判別すべき所望の画像に対応する画像データを読み出す。続くステップＳ３では、読み出した画像を小領域ａiに分割する。上述したように、小領域ａiは３２画素×３２画素の大きさであるため、分割される小領域ａiの個数は４０９６である。
【００４２】
続いて、ステップＳ５では、ＣＰＵはカウンタのカウント値を初期化（ｉ＝１）し、後で詳細に説明するように、ステップＳ７ではｉ番目の小領域ａiについてのテクスチャ特徴量Ｄ（第１特徴量），Ｖ（第２特徴量）およびＰ（第３特徴量）の計算処理を実行して、ステップＳ９でカウント値ｉが４０９６かどうかを判断する。つまり、すべての小領域ａiについてのテクスチャ特徴量を計算したかどうかを判断する。
【００４３】
ステップＳ９で“ＮＯ”であれば、つまりカウント値ｉが４０９６でなければ、すべての小領域ａiについてのテクスチャ特徴量を終了していないと判断して、ステップＳ１１でカウント値ｉをインクリメント（ｉ＝ｉ＋１）してからステップＳ７に戻る。一方、ステップＳ９で“ＹＥＳ”であれば、つまりカウント値ｉが４０９６であれば、すべての小領域ａiについてテクスチャ特徴量の計算処理を終了したと判断して、ステップＳ１３に進む。
【００４４】
ステップＳ１３では、後で詳細に説明するテクスチャ特徴量（Ｄ，Ｖ，Ｐ）を正規化し、正規化した特徴量（正規化特徴量）ＺＤ，ＺＶおよびＺＰを求める。続くステップＳ１５では、正規化特徴量ＺＤ，ＺＶおよびＺＰをＺＤ，ＺＶ，ＺＰ（３次元）空間上にプロットし、散布図（図５参照）を描く。
【００４５】
そして、ステップＳ１７で、正規化特徴量ＺＤ,ＺＶおよびＺＰの散布状態をＺＤ，ＺＶ，ＺＰ空間上で３個のクラスタにクラスタリングしてから処理を終了する。
【００４６】
図３は、図２に示したステップＳ７の小領域ａiのテクスチャ特徴量の計算処理を示すフロー図であり、特徴量の計算処理が開始されると、ＣＰＵは、ステップＳ２１で小領域ａiの画像データをＲ，Ｇ，ＢデータからＹデータに変換する。たとえば、小領域ａi上の位置（ｎ，ｍ）における画素データをｑi（ｎ，ｍ）とすると、各画素はＲ，Ｇ，Ｂの３種類の画素値を持つので、数８のように表記することができる。
【００４７】
【数８】
0003894802.tif 000009
【００４８】
ただし、ｎ，ｍは、それぞれ、１≦ｎ≦３２，１≦ｍ≦３２を満たす自然数である。ここで、位置すなわち座標系（ｎ，ｍ）は、画像に対する一般的な座標系（行列座標系）であり、水平位置ｎは右向きを正の方向とし、垂直位置ｍは下向きを正の方向とする。したがって、小領域ａiの画像の左上が（１，１）であり、右下が（３２，３２）である。
【００４９】
以後の処理においては、テクスチャは濃淡画像として扱うため、ＣＰＵは、数９に従って各画素のＲ，Ｇ，ＢデータをＹデータに変換する。
【００５０】
【数９】
0003894802.tif 000010
【００５１】
また、これ以降の処理では、テクスチャの平均的な明るさ（平均輝度）は問題にしないので、次のステップＳ２３では、Ｙデータから小領域ａiの平均輝度を除去する。
【００５２】
まず、数１０に従って小領域ａiの平均輝度を求める。
【００５３】
【数１０】
0003894802.tif 000011
【００５４】
次に、小領域ａiの各画素の輝度値Ｙi（ｎ，ｍ）から小領域ａiの平均輝度を減算した減算値（減算データ）ＹＴi（ｎ，ｍ）を数１１に従って求める。
【００５５】
【数１１】
0003894802.tif 000012
【００５６】
ただし、１≦ｎ≦３２,１≦ｍ≦３２
次に、ステップＳ２５では、小領域ａiの減算データＹＴi（ｎ，ｍ）を２次元離散的フーリエ変換（２次元ＤＦＴ）する。つまり、数１２に従って演算し、複素数Ｆi（ｕ，ｖ）を算出する。
【００５７】
【数１２】
0003894802.tif 000013
【００５８】
ただし、ｊは１単位虚数であり、ｕ，ｖは、それぞれ、−１６≦ｕ≦１６，−１６≦ｖ≦１６を満たす整数である。この複素数Ｆi（ｕ，ｖ）は３２画素×３２画素の減算データＹＴi（ｎ，ｍ）の２次元周波数（ｕ，ｖ）における周波数成分を表す。
【００５９】
ここで、ｕの単位は[cpw : cycle per picture width]であり、ｖの単位は[cph : cycle per picture hight]である。
一般的に、或る物理量が複素数で表される場合、その物理量の大きさや強さ（強度）は、その複素数の大きさすなわち絶対値で表すことができる。そこで、ＣＰＵは、続くステップＳ２７で２次元周波数平面上に数１２によって得られた複素数Ｆi（ｕ，ｖ）の絶対値を取る（算出する）。具体的には、数１３に従って複素数Ｆi（ｕ，ｖ）の絶対値（「Ｗi（ｕ，ｖ）」と表記する。）を取り、絶対値Ｗi（ｕ，ｖ）をもって２次元周波数成分Ｆi（ｕ，ｖ）の強度とする。
【００６０】
【数１３】
0003894802.tif 000014
【００６１】
ただし、上述したように、ｕ，ｖは、それぞれ、−１６≦ｕ≦１６，−１６≦ｖ≦１６を満たす整数である。
【００６２】
また、座標系（ｎ，ｍ）に対応して、２次元周波数座標系（ｕ，ｖ）についても、水平周波数ｕは右向きを正の方向とし、垂直周波数ｖは下向きを正の方向とする。このような２次元周波数座標系（平面上）に、２次元周波数強度分布Ｗi（ｕ，ｖ）を等高線表示した一例が図４のように示される。この図４では、グレースケールによって周波数分布の強度を表しており、黒に近いほど、強度が弱く、白に近いほど、強度が強い。
【００６３】
続いて、ステップＳ２９では、この強度分布Ｗi（ｕ，ｖ）を２次元周波数平面上の荷重分布とみなし、その重心位置を決定する（求める）。ただし、図４から分かるように、荷重分布（強度分布）は点（原点）対称であるため、重心位置が常に原点となることを避けるように、この実施例では、周波数平面の右半分の領域（水平周波数ｕが０から＋１６cpw、垂直周波数ｖが−１６から＋１６cphの範囲）についての重心Ｇiの位置（ｕiG，ｖiG）が数１４に従って求められる。
【００６４】
【数１４】
0003894802.tif 000015
【００６５】
ここで、ＳＷiは小領域ａiの周波数強度分布の右半分の領域における総荷重を表し、数１５で与えられる。また、重心位置の測定例が図４内に点Ｐで示される。
【００６６】
【数１５】
0003894802.tif 000016
【００６７】
なお、この実施例では、重心位置を周波数平面の右半分の領域で求めたが、これに限定する必要はなく、周波数平面の左半分、上半分或いは下半分の領域で求めてもよい。
【００６８】
また、計算の煩わしさを考慮しなければ、上、下、左、右半分の領域以外の場合で、原点を通る直線で周波数平面を２分割した一方の領域（範囲）内で重心位置を求めるようにしてもよい。
【００６９】
次に、この実施例では、テクスチャの方向は問題にしていないので、ステップＳ３１では、ＣＰＵは重心位置によらず，原点から重心位置までの距離Ｄiを小領域ａiのテクスチャの第１特徴量（テクスチャ特徴量Ｄ）として数１６により算出する。
【００７０】
【数１６】
0003894802.tif 000017
【００７１】
続いて、ステップＳ３３で、小領域aiのテクスチャの第２特徴量（テクスチャ特徴量Ｖ）として、上述したような周波数平面の右半分の領域の周波数範囲で、周波数強度分布の重心周りの分散ｖiを数１７により求める。このとき、周波数範囲の強度分布が、この小領域ａiに関する全サンプルであるので、不偏分散を求める必要はない。このため、数１７においては、分母がＳＷi-１ではなく、ＳＷiとされる。
【００７２】
【数１７】
0003894802.tif 000018
【００７３】
ここで、小領域のテクスチャに関して、右半分の領域における水平周波数ｕを第１変量、垂直周波数ｖを第２変量とすると、この２つの変量ｕ，ｖを次のように統合して新しい２つの変量、すなわち第１主成分zi1および第２主成分zi2に変換することができる。
【００７４】
まず、第１変量ｕの分散（水平周波数ｕに沿って分布する荷重による分散）Ｖiuu、第２変量ｖの分散（垂直周波数ｖに沿って分布する荷重による分散）Ｖivv、および第１変量ｕ、第２変量ｖの共分散Ｖiuvは数１８〜数２０によってそれぞれ求められる。
【００７５】
【数１８】
0003894802.tif 000019
【００７６】
【数１９】
0003894802.tif 000020
【００７７】
【数２０】
0003894802.tif 000021
【００７８】
ここで、これらの分散値を要素とする分散・共分散行列（実対称行列）Ｖiは数２１で示される。
【００７９】
【数２１】
0003894802.tif 000022
【００８０】
この実対称行列Ｖiの固有値λi1，λi2（ただし、λi1≧λi2≧０である。）および固有値λi1，λi2に対応する固有ベクトルｌi1，ｌi2は次のように求められる。具体的には、固有値λi1，λi2は、Ｖiの固有多項式（λに関する２次方程式）、すなわち数２２の根として求められる。また、固有ベクトルｌi1，ｌi2は、数２３を解くことにより求められる。
【００８１】
【数２２】
0003894802.tif 000023
【００８２】
【数２３】
0003894802.tif 000024
【００８３】
次に、新しい変量である第１主成分zi1および第２主成分zi2は、ｕ，ｖの１次結合であり、数２４によってそれぞれ定義される。
【００８４】
【数２４】
0003894802.tif 000025
【００８５】
ここで、固有ベクトルｌi1，ｌi2は、それぞれ、数２５のように表している。
【００８６】
【数２５】
0003894802.tif 000026
【００８７】
以上より、ＣＰＵは、ステップＳ３５で、２次元周波数平面の右半分の領域で第１主成分の寄与率Ｐiを小領域ａiのテクスチャの第３特徴量（テクスチャ特徴量Ｐ）として求める。すなわち、第１主成分の寄与率Ｐiは数２６で算出することができる。
【００８８】
【数２６】
0003894802.tif 000027
【００８９】
しかし、第１固有値λi1は、第１主成分zi1の分散Ｖimに等しく、第２固有値λi2は、第２主成分zi2の分散Ｖisに等しいことが知られており、また、この実施例では、全変量の個数は２個であることから、第１主成分ｚi1の分散Ｖimと第２主成分ｚi2の分散Ｖisの和は、全分散Ｖiすなわち重心周りの分散に等しくなる。このため、寄与率Ｐiは、第１主成分ｚi1の分散Ｖimが全分散Ｖiに占める割合であるということもでき、数２７が成立する。
【００９０】
【数２７】
0003894802.tif 000028
【００９１】
なお、この第１主成分ｚi1の座標軸は、重心を通るあらゆる軸の中で、その軸に沿った分散が最も大きくなる軸を意味する。第１主成分軸と第２主成分軸との測定例は、図４においてそれぞれ実線と点線とで示される。
【００９２】
以上のように、Ｓ３１，Ｓ３３およびＳ３５の処理を経て、この小領域ａiのテクスチャ特徴量すなわち、距離Ｄi、分散Ｖi、第１主成分の寄与率Ｐiが得られる。
【００９３】
このような３種類のテクスチャ特徴量Ｄi，Ｖi，Ｐiは、それぞれについて４０９６個の測定値が得られる。したがって、上述したステップＳ１３においては、これら３種類の特徴量Ｄi，Ｖi，Ｐiのそれぞれについて、その４０９６個の測定データに渡る平均値が０、分散が１になるように正規化が行なわれる。以下では、正規化された特徴量（正規化特徴量）を、それぞれ、ＺＤi，ＺＶi，ＺＰiと表記することにする。
【００９４】
まず、テクスチャ特徴量Ｄiは、数２８に従って正規化される。
【００９５】
【数２８】
0003894802.tif 000029
【００９６】
ここで、ＤM，σDは、それぞれＤiの平均および標準偏差を表し、数２９および数３０によって計算される。
【００９７】
【数２９】
0003894802.tif 000030
【００９８】
【数３０】
0003894802.tif 000031
【００９９】
また、テクスチャ特徴量Ｖiは、数３１に従って正規化される。
【０１００】
【数３１】
0003894802.tif 000032
【０１０１】
ここで、ＶM，σVは、それぞれＶiの平均および標準偏差を表し、数３２および数３３によって計算される。
【０１０２】
【数３２】
0003894802.tif 000033
【０１０３】
【数３３】
0003894802.tif 000034
【０１０４】
さらに、テクスチャ特徴量Ｐiは、数３４に従って正規化される。
【０１０５】
【数３４】
0003894802.tif 000035
【０１０６】
ここで、ＰM，σPは、それぞれＰiの平均および標準偏差を表し、数３５および数３６によって計算される。
【０１０７】
【数３５】
0003894802.tif 000036
【０１０８】
【数３６】
0003894802.tif 000037
【０１０９】
以下、この実施例において、４０９６個の小領域ａiのテクスチャ特徴量Ｄi，Ｖi，Ｐiをまとめて（Ｄ，Ｖ，Ｐ）と表記することとする。同様に、正規化特徴量ＺＤi，ＺＶi，ＺＰiをまとめて（ＺＤ，ＺＶ，ＺＰ）と表記することにする。
【０１１０】
また、上述したステップＳ１５では、４０９６個の正規化特徴量の測定値（ＺＤ，ＺＶ，ＺＰ）が、ＺＤ，ＺＶ，ＺＰ空間上にプロットされる。ただし、この実施例では、簡単に説明するため、或る画像について、１８個の測定値（ＺＤ，ＺＶ，ＺＰ）がプロットされた例を図５に示す。
【０１１１】
本件出願人の実験結果によれば、模様が平坦なテクスチャであるほど、距離Ｄ、分散Ｖ、第１主成分の寄与率Ｐが大きく、したがって、ＺＤ，ＺＶ，ＺＰが大きい。逆に、模様が明確なテクスチャであるほど、距離Ｄ，分散Ｖ，第１主成分の寄与率Ｐが小さく、したがって、ＺＤ，ＺＶ，ＺＰが小さい。このように、著しい特徴があることが判っている。
【０１１２】
このことは、理論的にも説明することができる。まず、模様が平坦なテクスチャであるほど距離Ｄが大きく、模様が明確なテクスチャであるほど距離Ｄが小さくなるのは、▲１▼平坦なテクスチャほど（輝度の平均値を除算した後では）、主たる強度の集中がなく、同程度の小さな強度のみが広がるので、強度分布が完全に一様な場合の重心位置Ｇ０（ほぼ水平周波数が８cpw，垂直周波数が０cph）に近づくこと、▲２▼この実施例で扱っているような自然画像では、強度の集中個所はかなり原点に近いこと、したがって、非常に平坦なテクスチャであっても、その重心位置は図中のＧ０で示す位置よりも原点寄りにあること、▲３▼また、模様が明確であるほど、重心位置はより強く強度の集中個所に引きずられて、図中のＧ０で示す位置からより多く原点寄りになることによるためである。
【０１１３】
このことから、テクスチャが平坦であるかどうかは第１特徴量Ｄに大きく依存していると言える。したがって、少なくとも第１特徴量Ｄを取得することができれば、平坦なテクスチャかどうかを判別することができる。
【０１１４】
また、第１特徴量Ｄが小さいと第２特徴量Ｖも小さくなる理由は次のとおりである。模様が明確であればあるほど、Ｄは小さくなるが、その重心付近に大きな強度の集中が現れるので、模様が平坦な場合に比べて，まず第１主成分の分散Ｖmが大幅に減少する。つまり、分散は、重心から離れたところに主たる重量があると急激に大きくなり、逆に重心近くにのみ主たる重量があると急激に小さくなる。また、第２主成分の分散Ｖsも、分散Ｖmほどは大幅には減少しないが、減少する。分散Ｖは、分散Ｖmと分散Ｖsとの和であるから、結果として、模様が明確であればあるほど、第１特徴量Ｄおよび第２特徴量Ｖが小さくなる。逆に模様が平坦であればあるほど、第１特徴量Ｄおよび第２特徴量Ｖが大きくなる。
【０１１５】
さらに、第２特徴量Ｖが小さいと第３特徴量Ｐも小さくなる理由ついては次のとおりである。模様が明確であれば、上述のとおり、分散Ｖmとともに分散Ｖsも減少するが、分散Ｖsの減少の度合いは、分散Ｖmの減少の度合いよりも緩やかである。したがって、模様が明確であるほど、分散Ｖsの分散Ｖmに対する比率が相対に増えることになる。すなわち、分散Ｖsを数３７のように表すと、模様が明確であるほど、αが大きくなる。したがって、第３特徴量Ｐは、数３８に示すように、小さくなる。
【０１１６】
【数３７】
Ｖs=α・Ｖm
【０１１７】
【数３８】
Ｐ＝Ｖm／（Ｖm＋Ｖs）＝１／（１＋α）
結果として、模様が明確であればあるほど、第２特徴量Ｖおよび第３特徴量Ｐが小さくなる。逆に模様が平坦であればあるほど、第２特徴量Ｖも第３特徴量Ｐも大きくなる。
【０１１８】
以上のように、模様が平坦なテクスチャであるほどテクスチャ特徴量（Ｄ，Ｖ，Ｐ）が大きく、模様が明確なテクスチャであるほどテクスチャ特徴量（Ｄ，Ｖ，Ｐ）が小さいということが理論的に説明できる。
【０１１９】
したがって、この実施例では、上述したステップＳ１７で説明したように、ステップＳ１５においてＺＤ，ＺＶ，ＺＰ空間上にプロットされた分布図に対し、ＺＤ，ＺＶ，ＺＰが共に最も小さいクラスタおよびＺＤ，ＺＶ，ＺＰが共に最も大きいクラスタをそれぞれ一まとめにしてクラスタリングを行ない、残りの点をまとめて１クラスタとするようにしてある。このようにして、画像内の４０９６個の小領域ａiが、模様のかなり平坦なもの、模様のかなり明確なもの、その中間的なものの３つのクラスタに分類される。
【０１２０】
このとき、一般に流通している統計処理ソフトウェアを用いて、クラスタ間の距離の決め方としてたとえば重心法によるクラスタリングを実行すれば、容易にクラスタリングできる。
【０１２１】
なお、一般的に流通している統計処理ソフトウェアとしては、ＳＰＳＳ  ｌｎｃ．社製の「ＳＰＳＳ」やＳｔａｔＳｏｆｔ社製の「ＳＴＡＴＩＳＴＩＣＡ」を用いることができる。
【０１２２】
この実施例によれば、テクスチャ特徴量Ｄ、ＶおよびＰという３種類の物理量の組み合わせを用いるようにし、さらにその組み合わせが画像全体或いは対象とする集合全体に渡るテクスチャ特徴量Ｄ、ＶおよびＰの分布状況のどこに位置するかによりテクスチャ平坦度を判別するようにしたので、判別精度を向上させることができる。
【０１２３】
また、小領域（３２画素×３２画素）の全部に渡って抽出された物理量を求めるようにしたので、画素単位に重畳しているノイズの影響は、そのランダム的性質により、３２画素×３２画素に渡る演算の過程で相殺され、精度の良い物理量が得られる。
【０１２４】
なお、上述の実施例では、小領域ａiに対する前処理としてＲ，Ｇ，ＢデータからＹデータを求めるようにしたが、Ｒ，Ｇ，Ｂデータから色相Ｈ，彩度Ｓ，明度Ｉ（または、Ｖと表記されることもある。）のデータすなわちＨＳＩ空間のデータに変換し、その結果である明度Ｉのデータを用いるようにしてもよい。
【０１２５】
また、上述の実施例では、小領域ａiを３２画素×３２画素の大きさとしたが、必ずしもこの大きさに限定する必要はなく、対象とするテクスチャとして最も分析しやすい大きさを任意に選べばよい。
【０１２６】
さらに、上述の実施例では、２次元周波数成分の算出にあたり、２次元離散的フーリエ変換（ＤＦＴ）を用いたが、特に２次元高速フーリエ変換（ＦＦＴ）、或いは、２次元離散的コサイン変換（ＤＣＴ）を用いてもよい。
【０１２７】
さらにまた、上述の実施例では、３個のクラスタにクラスタリングする場合を示したが、これに限る必要はなく、目的に応じて２個或いは５個というように増減してもよいことはもちろんである。
【０１２８】
次に本発明の第２の実施例について説明する。この第２の実施例のテクスチャ判別システムは、テクスチャ特徴量のクラスリング方法が異なる以外は上述の実施例と同じであるため、重複した説明は省略する。
【０１２９】
つまり、クラスタ数を最適なクラスタ数（最適クラスタ数）に決定し、その最適クラスタ数でクラスタリングを実行するようにしてある。ここで、クラスタリング処理手順を簡単に説明することにする。まず、クラスタ数ａがＮ個（全標本個数と同じ数）の場合に、クラスタ間の距離の決め方としてたとえば重心法に基づいてクラスタリングすると、クラスタリングは一意に決まる。
【０１３０】
クラスタ数がａ（ただし、ａ＝Ｎ，Ｎ−１，…，２，１）のときのｉ番目のクラスタをＣ（ａ，ｉ）と表記する。ここで、ｉ＝１，２，…，ａである。また、クラスタＣ（ａ,ｉ）に含まれる標本点の個数をｂ（ａ，ｉ）とし、クラスタＣ（ａ，ｉ）に含まれるｋ番目の標本点をＳ（ａ，ｉ，ｋ）と表記する。ここで、ｋ＝１，２，…，ｂ（ａ，ｉ）である。さらに、クラスタＣ（ａ，ｉ）の重心をＧ（ａ，ｉ）とし、重心Ｇ（ａ，ｉ）と標本点Ｓ（ａ，ｉ，ｋ）との距離をｄ（ａ，ｉ，ｋ）と表記すると、クラスタＣ（ａ，ｉ）内の各標本点を重心Ｇ（ａ，ｉ）で代表させたことによるクラスタＣ（ａ，ｉ）についてのクラスタリング誤差は、数３９で示される。
【０１３１】
【数３９】
0003894802.tif 000038
【０１３２】
なお、以下においては、簡単のため、このＤ（ａ，ｉ）をｉ番目のクラスタＣ（ａ，ｉ）の誤差と呼ぶことにする。
【０１３３】
最初のクラスタ数ａ＝Ｎのときには、各クラスタＣ（Ｎ，ｉ）（ただし、ｉ＝１,２,…,Ｎ）は1個の標本点のみを含むので、各クラスタＣ（Ｎ，ｉ）内の標本点と各クラスタの重心Ｇ（Ｎ，ｉ）は一致する。したがって、ｂ（Ｎ，ｉ）＝１，ｋ＝１，ｄ（Ｎ，ｉ，ｋ）＝０となるので、Ｄ（ａ，ｉ）＝０となる。すなわち、ａ＝Ｎのとき、どのクラスタＣ（ａ，ｉ）についても、クラスタＣ（ａ，ｉ）の誤差は０である。
【０１３４】
また、ｉ番目のクラスタＣ（ａ，ｉ）の誤差Ｄ（ａ，ｉ）を（クラスタ数がａのときの）全クラスタ（ｉ＝１，２，…，ａ）について加算した誤差Ｄ（ａ，ｉ）の総和ＳＤ（ａ）は、数４０で算出される。
【０１３５】
【数４０】
0003894802.tif 000039
【０１３６】
つまり、ＳＤ（ａ）は、Ｎ個の標本点をａ個のクラスタにクラスタリングしたことによるクラスタリングの誤差を表していると考えられる。以下、第２の実施例では、このＳＤ（ａ）をａ個へのクラスタリングの誤差と呼ぶことにする。最初のクラスタ数ａ＝Ｎのときのクラスタリングについては、上述したように各クラスタの誤差Ｄ（ａ,ｉ）は０であるので、ＳＤ（ａ）＝ＳＤ（Ｎ）＝０となる。すなわちクラスタ数ａ＝Ｎのとき、ａ個へのクラスタリングの誤差は０である。
【０１３７】
次に、クラスタ数ａ＝Ｎ−１の場合に、重心法に基づいてＮ個の全標本点をクラスタリングすると、Ｎ個の標本点の内、最も近い２点をまとめて１つのクラスタとし、残りのクラスタはクラスタ数ａ＝Ｎの場合のままとして、クラスタ数ａ＝Ｎ−１個へのクラスタリングを決定する。そして、決定されたクラスタリングについて、上記の諸量を求める処理（数４０を用いた誤差の総和の算出処理）を繰り返せば、最も誤差の少ないクラスタ数ａ＝Ｎ−１個へのクラスタリングの誤差ＳＤ（ａ）＝ＳＤ（Ｎ−１）を得ることができる。
【０１３８】
続いて、クラスタ数ａ＝Ｎ−２の場合には、重心法に基づいて、現在のクラスタ数ａ＝Ｎ−１個のクラスタリングがクラスタ数ａ＝Ｎ−２個へのクラスタリングに変更される。この場合、Ｎ−１個の重心の内、最も近い２つの重心を１つのクラスタにまとめるように、クラスタ数ａ＝Ｎ−２個へのクラスタリングを決定する。
【０１３９】
なお、このようなクラスタリングの分け方の列挙は、一般に流通している統計解析ソフトウェアを使用するか、自分でプログラムを作成すれば、可能であり、さらに自動化することも容易である。以下、クラスタリングする場合については同様のことが言える。
【０１４０】
そして、決定されたクラスタリングで、上記の諸量を求める処理（数４０を用いた誤差の総和の算出処理）を繰り返せば、実用上差し支えのない程度に誤差の少ないクラスタ数ａ＝Ｎ−２個へのクラスタリングの誤差ＳＤ（ａ）＝ＳＤ（Ｎ−２）を得ることができる。
【０１４１】
ただし、Ｎ−２個へのクラスタリングを決める場合には、前回のクラスタリングの結果、すなわちＮ−１個へのクラスタリングの結果を用いずに、Ｎ個の標本をＮ−２個のクラスタに分けるすべての分け方について、上記の手順によってＳＤ（ａ）＝ＳＤ（Ｎ−２）を算出し、その中で最も小さいものを真のＳＤ（ａ）とするようにしてもよい。
【０１４２】
この後者の方法の方がＳＤ（ａ）の評価としてはより正確ではあるが、全標本の個数Ｎが非常に大きい場合には、コンピュータの処理能力にもよるが、現実的な方法ではなくなってしまう。このような場合には、前者の方法を用いて、実用上差し支えのない程度に誤差の少ないＳＤ（ａ）を求めるのがよいと考えられる。
【０１４３】
ここで、クラスタに分けるすべての分け方について具体的に説明することにする。たとえば、５個（Ｎ＝５）の標本点を３個のクラスタに分ける（クラスタリングする）すべての分け方は次のように、場合分けすることができる。ただし、５個の標本点（要素）は１〜５の数字で表すことにする。
【０１４４】
（ｉ）要素の数が３，１，１の３つのクラスタに分ける場合には、クラスタの組み合わせは次のようになる。なお、組み合わせ数は１０（5Ｃ3）である。
（１，２，３），（４），（５）
（１，２，４），（３），（５）
（１，２，５），（３），（４）
（１，３，４），（２），（５）
（１，３，５），（２），（４）
（１，４，５），（２），（３）
（２，３，４），（１），（５）
（２，３，５），（１），（４）
（２，４，５），（１），（３）
（３，４，５），（１），（２）
（ii）要素の数が２，２，１の３つのクラスタに分ける場合には、クラスタの組み合わせは次のようになる。
（１，２），（３，４），（５）
（１，３），（２，４），（５）
（１，４），（２，３），（５）
（１，２），（３，５），（４）
（１，３），（２，５），（４）
（１，５），（２，３）， (４)
（１，２），（４，５），（３）
（１，４），（２，５），（３）
（１，５），（２，４），（３）
（１，３），（４，５），（２）
（１，４），（３，５），（２）
（１，５），（３，４），（２）
（２，３），（４，５），（１）
（２，４），（３，５），（１）
（２，５），（３，４），（１）
以上、（ｉ），（ii）で示すように、５個の標本点を３つのクラスタにクラスタリングする場合には、２５通りのクラスタリング方法がある。このように、分け方は多数存在するため、上述したように、前回の結果を利用してクラスタリングすることにより、コンピュータ（ＣＰＵ）の処理量を軽減してもよい。
【０１４５】
さらに、クラスタ数ａ＝Ｎ−３の場合には、重心法に基づいて、現在のクラスタ数ａ＝Ｎ−２個のクラスタリングがクラスタ数ａ＝Ｎ−３個へのクラスタリングに変更される。この場合、Ｎ−２個の重心の内、最も近い２つの重心を１つのクラスタにまとめるように、クラスタ数ａ＝Ｎ−３個へのクラスタリングを決め、上記の諸量を求める処理（数４０を用いた誤差の総和の算出処理）を繰り返せば、実用上差し支えのない程度に誤差の少ないクラスタ数ａ＝Ｎ−３個へのクラスタリングの誤差ＳＤ（ａ）＝ＳＤ（Ｎ−３）を得ることができる。
【０１４６】
ただし、上述した場合と同様に、Ｎ−３個へのクラスタリングを決める場合には、前回のクラスタリングの結果、すなわちＮ−２個へのクラスタリングの結果を用いずに、Ｎ個の標本をＮ−３個のクラスタに分けるすべての分け方について上記の手順にてＳＤ（ａ）＝ＳＤ（Ｎ−３）を算出し、その中で最も小さいものを真のＳＤ（ａ）とするようにしてもよい。
【０１４７】
このようにして、クラスタリングが繰り返され、クラスタ数ａ＝Ｎ−（Ｎ−２）＝２の場合には、重心法に基づいて、現在のクラスタ数ａ＝３個のクラスタリングが２個のクラスタリングに変更される。この場合、３個の重心の内、最も近い２つの重心を１つのクラスタにまとめるように、クラスタ数ａ＝２個へのクラスタリングを決め、上記の諸量を求める処理（数４０を用いた誤差の総和の算出処理）を繰り返せば、実用上差し支えのない程度に誤差の少ないクラスタ数（ａ＝２）へのクラスタリングの誤差ＳＤ（ａ）＝ＳＤ（２）を得ることができる。
【０１４８】
ただし、２個へのクラスタリングを決める場合には、前回のクラスタリングの結果、つまり、３個へのクラスタリングの結果を用いずに、Ｎ個の標本を２個のクラスタに分けるすべての分け方について、上記の手順によってＳＤ（ａ）＝ＳＤ（２）を算出し、その中で最も小さいものを真のＳＤ（ａ）とするようにしてもよい。
【０１４９】
最後に、ａ＝Ｎ−（Ｎ−１）＝１の場合には、重心法に基づいて、現在のクラスタ数ａ＝２個のクラスタリングをクラスタ数ａ＝１個へのクラスタリングに変える。この場合のクラスタリングは一意に決まる。統合された1つのクラスタの重心は２個の重心の中点となる。このクラスタリングについて、上記の諸量を求める処理（数４０を用いた誤差の総和の算出処理）を繰り返せば、最も誤差の少ないａ＝１個へのクラスタリングの誤差ＳＤ（ａ）＝ＳＤ（１）を得ることができる。
【０１５０】
上述のようにして得たａ個へのクラスタリングの誤差ＳＤ（ａ）（ただし、ａ＝Ｎ，Ｎ−１，…，２，１）を、横軸をクラスタ数ａとしてプロットする。第２の実施例では、Ｎ＝１８であり、前回のクラスタリングの結果を用いる方法で得たＳＤ（ａ）をプロットした一例を図６に示す。
【０１５１】
図６を参照して、最適なクラスタ数（最適クラスタ数）ａoptは、クラスタ数ａをＮから１に向かって減少させたときに、急激にＳＤ（ａ）が増大（変化）する個所を含む変化前（直前）のクラスタ数ａ、または、それより大きなクラスタ数ａであって、許せる程度の大きさ（許容範囲）のＳＤ（ａ）を与えるクラスタ数ａに決定する。したがって、図５に示す例では、最適クラスタ数ａoptは「３」である。
【０１５２】
ただし、この実施例では、許容範囲は最大誤差の２０パーセント未満としてあるが、設計者や使用者等によって自由に変更可能である。また、最適クラスタ数ａoptは、急激にＳＤ（ａ）が変化する変化後（直後）のクラスタ数ａであって、許容範囲を超えないクラスタ数ａに決定するようにしてもよい。したがって、許容範囲の決め方によっては、数個所の最適クラスタ数ａoptの候補が存在することとなり、その中から１つを決定する場合がある。
【０１５３】
上述したようなクラスタリング処理は、コンピュータのＣＰＵが図７および図８に示すフロー図に従って処理する。なお、第２の実施例におけるテクスチャ平坦度判別処理では、図２に示したステップＳ１３〜Ｓ１７の処理に代えて、図７および図８に示す処理が実行される。
【０１５４】
具体的には、クラスタリング処理を開始すると、ＣＰＵはステップＳ４１でクラスタ数ａを初期化（ａ＝Ｎ）する。続いて、ステップＳ４３では、重心法によりＮ個の全標本値をクラスタリングする。そして、ステップＳ４５でカウンタのカウント値ｉを初期化（ｉ＝１）して、ステップＳ４７でクラスタＣ（ａ，ｉ）の重心Ｇ（ａ，ｉ）を算出する。
【０１５５】
ステップＳ４９では、カウンタのカウント値ｋを初期化（ｋ＝１）し、ステップＳ５１では、クラスタＣ（ａ，ｉ）内の重心Ｇ（ａ，ｉ）と各標本点Ｓ（ａ，ｉ，ｋ）との距離ｄ（ａ，ｉ，ｋ）を算出する。
【０１５６】
次に図８に示すように、次のステップＳ５３では、ｋ＝ｂ（ａ，ｉ）かどうかを判断する。つまり、クラスタＣ（ａ，ｉ）における全標本点との距離ｄ（ａ，ｉ，ｋ）を算出したかどうかを判断する。
【０１５７】
ステップＳ５３で“ＮＯ”であれば、つまりｋ＝ｂ（ａ，ｉ）でなければ、クラスタＣ（ａ，ｉ）における全標本点との距離ｄ（ａ，ｉ，ｋ）を算出していないと判断し、ステップＳ５５でカウント値ｋをインクリメント（ｋ＝ｋ＋１）して、図７に示したステップＳ５１に戻る。
【０１５８】
一方、ステップＳ５３で“ＹＥＳ”であれば、つまりｋ＝ｂ（ａ，ｉ）であれば、クラスタＣ（ａ，ｉ）における全標本点との距離ｄ（ａ，ｉ，ｋ）を算出した判断し、ステップＳ５７でｉ番目のクラスタＣ（ａ，ｉ）の誤差Ｄ（ａ，ｉ）を数３９に従って算出する。
【０１５９】
続くステップＳ５９では、カウント値ｉがクラスタ数ａと等しいかどうかを判断する。このステップＳ５９で“ＮＯ”であれば、つまりカウント値ｉがクラスタ数ａと等しくなければ、ステップＳ６１でカウンタをインクリメント（ｉ＝ｉ＋１）して、図７に示したステップＳ４７に戻る。
【０１６０】
一方、ステップＳ５９で“ＹＥＳ”であれば、つまりカウント値ｉがクラスタ数ａと等しければ、ステップＳ６３でａ個へのクラスタリングの誤差ＳＤ（ａ）を算出する。具体的には、数４０を用いて、ａ個のクラスタにクラスタリングしたことによるクラスタリングの誤差を求める。
【０１６１】
続いて、ステップＳ６５では、クラスタ数ａが１かどうかを判断する。このステップＳ６５で“ＮＯ”であれば、つまりクラスタ数ａが１でなければ、ステップＳ６７でクラスタ数ａをデクリメントし（ａ＝ａ−１）し、重心法に基づいて、今回の（デクリメント後の）クラスタ数ａ（ａ個）のクラスタにクラスタリングしてから、図７に示したステップＳ４５に戻る。
【０１６２】
一方、ステップＳ６５で“ＹＥＳ”であれば、つまりクラスタ数ａが１であれば、ステップＳ６９で誤差ＳＤ(ａ)が急激に増大する直前のクラスタ数ａ、または、それより大きなクラスタ数ａで許容できる大きさ（許容範囲）における誤差ＳＤ（ａ）を与えるａを最適クラスタ数ａoptに決定してから処理を終了する。
【０１６３】
ただし、第２の実施例では、誤差ＳＤ(ａ)の許容範囲は、最大誤差（約５５）の１／５（約１０）とし、誤差ＳＤ（ａ）が１０を超えない範囲で最適クラスタ数ａoptが決定される。つまり、許容範囲は、プログラム（ソフト）の設計者や使用者が任意に設定できる値である。
【０１６４】
第２の実施例によれば、クラスタ数を妥当な値に決定してクラスタリングできるので、単なる統計として検証する際の利便性を向上することができるとともに、クラスタリングによってその後の処理に悪影響を与えるのを防止することができる。また、このように、最適なクラスタ数にクラスタリングするようにすれば、上述の実施例のように、明確な模様、平坦な模様およびその中間のような単純な分け方でなく、さらに複数段階の分け方が可能である。
【０１６５】
なお、第２の実施例では、クラスタ数ａを１個ずつ減らすようにしたが、全標本個数Ｎが膨大である場合には、クラスタ数ａを複数個（たとえば、１０個）ずつ減らすようにしてもよい。
【０１６６】
また、第２の実施例では、クラスタリングの手法として、ユークリッド距離に基づく重心法を適用した場合についてのみ説明したが、類似度を表す距離としてはユークリッド距離に限定する必要はない。たとえば、ユークリッド距離の２乗、マハラノビス距離、相関係数等を用いるようにしてもよい。
【０１６７】
さらに、クラスタ間の距離の決め方についても、重心法に限定する必要はなく、最短距離法、最長距離法、群平均法、メディアン法またはウォード法等を用いることも可能である。
【図面の簡単な説明】
【図１】この発明の一実施例のテクスチャ平坦度判別システムで取り扱われる画像を示す図解図である。
【図２】テクスチャ平坦度判別を示すフロー図である。
【図３】テクスチャ特徴量の計算処理を示すフロー図である。
【図４】テクスチャ特徴量の計算処理において、減算データを２次元ＤＦＴして得られた複素数の絶対値を２次元周波数平面上に等高線表示した図解図である。
【図５】テクスチャ平坦度判別処理において、テクスチャ特徴量に従ってプロットしたグラフである。
【図６】この発明の第２の実施例におけるクラスタリング処理のクラスタ数に対する誤差を示すグラフである。
【図７】この発明の第２の実施例のクラスタリング処理の一部を示すフロー図である。
【図８】この発明の第２の実施例のクラスタリング処理の他の一部を示すフロー図である。
【図１】
0003894802.tif 000040
【図２】
0003894802.tif 000041
【図３】
0003894802.tif 000042
【図４】
0003894802.tif 000043
【図５】
0003894802.tif 000044
【図６】
0003894802.tif 000045
【図７】
0003894802.tif 000046
【図８】
0003894802.tif 000047