x-means

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.19.3377

クラスタリングの方法は、色々あるけども、わたしの知る限り、ほとんどの方法では、事前にいくつのクラスターに分かれるか決める必要がある。クラスタリングしたい多くの状況で、いくつのクラスターに分かれるか事前には分からないというシチュエーションは山ほどあると思う。

ネットワークのクラスタリングには、Newman法とかあるらしいけど、x-meansはk-meansの最適なkを自動で決定するアルゴリズムであるらしい。やってることは単純で、k=2でk-meansを使って、どんどん分割していって、適切なとこまで分割できたら、そこで分割を打ち切る。こんなことは、最適なクラスター数が分からんという事態に直面したら誰でも考えるので、適切かどうかの判定をどうやってやるかが肝。論文ではBICないしAICというのを使えばいいよと書いてある。BICが何かは知らないけれど、対数尤度に謎の項を足したもの。多分大雑把には、k-meansは各クラスターがクラスター中心の周りに正規分布していると仮定しているようなものなので、このモデルがどれくらい適合しているかを測る指標があればよいという感じなんじゃないかと思う。まあ、BICが何か分からなくても、論文に具体的に式が書いてあるので実装はできる

OpenCVのサンプルに、k-meansを使った減色/画像分割があるので、これをx-meansで置き換えてみる
http://opencv.jp/sample/misc.html#color-sub
実際にやってみると、クラスター数が数百になったりして(画像サイズにもよるだろうけど)、相当遅い。ていうか、実装は合ってるのか怪しい。まあ、結論としては、意味なかった。思うに、共分散が単位行列に比例するとしているので、特定方向のみ分散が大きいようなケースでは分割しすぎるという問題があるんじゃなかろうか。あと、情報量基準でよいモデルを選択できるという根拠がいまいちよくわかってない。もっとroughに適合度検定とかじゃダメなんだろうか

実装。OpenCVの新しいPythonインターフェースは分からないことが多いので、とりあえず動くことだけを考えて作ってる

import sys
sys.path.append("C:\\\\OpenCV2.1\\Python2.6\\Lib\\site-packages")

from cv import *
import numpy
from numpy.linalg import norm
from math import log
from struct import *

def xmeans(vectors , N0=1):
    R = len(vectors)
    M = len(vectors[0])
    def IC(indices , K):
       eps = 1.0e-6
       Rn = float(len(indices))
       center= reduce(lambda x,y:x+y , [vectors[idx] for idx in indices] , numpy.zeros(M))/Rn
       variance = eps + sum([norm(vectors[idx]-center)**2 for idx in indices])/Rn
       #AIC
       #ic = Rn*log(2*numpy.pi) + Rn*M*log(variance) + (Rn-K) - 2*Rn*log(Rn/R) + 2*(M+1)*K
       #BIC
       ic = Rn*log(2*numpy.pi) + Rn*M*log(variance) + (Rn-K) - 2*Rn*log(Rn/R) + K*(M+1)*log(Rn)
       return ic
    def kmeans(indices , N):
       clusters = CreateMat(len(indices) , 1 , CV_32SC1)
       points = CreateMat(len(indices) , M  , CV_32FC1)
       for i in xrange(points.rows):
          for j in xrange(points.cols):
             idx = indices[i]
             points[i,j] = vectors[idx][j]
       KMeans2(points , N ,  clusters , (CV_TERMCRIT_EPS + CV_TERMCRIT_ITER, 10, 1.0))
       ret = [[] for _ in xrange(N)]
       for i in xrange(clusters.rows):
          idx = int(clusters[i,0])
          ret[idx].append(indices[i])
       return ret
    def split(indices):
       result = []
       stack = [indices]
       while True:
          if len(stack)==0:return result
          targets = stack[-1]
          if len(targets)<=1:
             stack.remove(targets)
             result.append(targets)
             continue
          ic1 = IC(targets , 1)
          D1,D2 = kmeans(targets,2)
          if len(D1)==0 or len(D2)==0:
             stack.remove(targets)
             result.append(targets)
             continue
          ic2 = IC(D1,2)+IC(D2,2)
          if ic2 < ic1:
             stack.remove(targets)
             stack.append(D1)
             stack.append(D2)
          else:
             stack.remove(targets)
             result.append(targets)
    result = []
    ls = kmeans(range(R) , N0)
    for indices in ls:
       result += split(indices)
    ret = numpy.zeros(R , numpy.uint32)
    for i,indices in enumerate(result):
       for idx in indices:
           ret[idx] = i
    return ret,len(result)

if __name__=="__main__":
  src_img = LoadImage("test.jpg")
  sz = src_img.width*src_img.height
  dst_img = CloneImage(src_img)
  vectors = []
  imgstr = src_img.tostring()
  for i in xrange(sz):
    c1,c2,c3 = unpack('BBB' , imgstr[i*3:i*3+3])
    vectors.append( numpy.array([float(c1),float(c2),float(c3)]) )

  print "start x-means"
  clusters,k = xmeans(vectors,2)
  print "%d clusters" % k

  count = numpy.zeros(k)
  r = numpy.zeros(k)
  g = numpy.zeros(k)
  b = numpy.zeros(k)
  for i in xrange(sz):
    idx = clusters[i]
    count[idx] += 1
    (cb,cg,cr) = tuple(vectors[i])
    r[idx] += cr
    g[idx] += cg
    b[idx] += cb
  for idx in range(k):
    if count[idx]>0.0:
      r[idx] /= count[idx]
      g[idx] /= count[idx]
      b[idx] /= count[idx]
  for h in xrange(dst_img.height):
    for w in xrange(dst_img.width):
      idx = clusters[h*dst_img.width+w]
      dst_img[h,w] = (b[idx] , g[idx] , r[idx])
  NamedWindow ("result", CV_WINDOW_AUTOSIZE)
  ShowImage ("result", dst_img)
  WaitKey (0)
  DestroyWindow("result")