まえおき

先日投稿した「実践コンピュータビジョン」の全演習問題をやってみたの詳細です。なお、あくまで「やってみた」であって、模範解答ではありません。「やってみた」だけなので、うまくいかなかった、という結果もあります。

この記事は「実践コンピュータビジョン」（オライリージャパン、Jan Erik Solem著、相川愛三訳）の演習問題に関する記事です。この記事を書いた人間と書籍との間には何の関係もありません。著者のJan氏と訳者の相川氏にはコンピュータビジョンについて学ぶ機会を与えて頂いたことに感謝します。

実践コンピュータビジョンの原書("Programming Computer Vision with Python")の校正前の原稿は著者のJan Erik Solem氏がCreative Commons Licenseの元で公開されています。。

回答中で使用する画像の多く、また、sfm.py, sift.py, stereo.pyなどの教科書中で作られるファイルはこれらのページからダウンロードできます。
オライリージャパンのサポートページ
 Jan Erik Solem氏のホームページ
 相川愛三氏のホームページ

ここで扱うプログラムはすべてPython2.7上で動作を確認しています。一部の例外を除きほとんどがJupyter上で開発されました。また必要なライブラリーを適宜インストールしておく必要があります。Jupyter以外の点は、上記の書籍で説明されています。

結構な分量なので章毎に記事をアップしています。今回は第６章です。

６章　画像のクラスタリング

いろいろなクラスタリング法を学びます。応用が効きそうです

6.4.1 Hierarchical k-means is a clustering method that applies k-means recursively to the clusters to create a tree of incrementally refined clusters. In this case, each node in the tree will branch to k child nodes. Implement this and try it on the font images.

階層型k平均法（k平均法を再帰的に適用したもの）をフォント画像に適用せよ

私の回答

ようするにk平均法を何回も適用して、それぞれの集合がk以下になったら止めればいいわけですね？たぶん。
imtools.pyやpickle.pyなど教科書内で使用するライブラリーが必要です¹。
また、フォントとその主成分分析結果を保存したもの（第１章で作成した"font_pca_modes.pkl"）が必要です。

from PIL import Image, ImageDraw
from numpy import *
from pylab import *
import scipy.misc
from scipy.cluster.vq import *
import imtools
import pickle

imlist = imtools.get_imlist('selected_fontimages/')
imnbr = len(imlist)
with open('font_pca_modes.pkl', 'rb') as f:
    immean = pickle.load(f)
    V = pickle.load(f)
immatrix = array([array(Image.open(im)).flatten() for im in imlist], 'f')
immean = immean.flatten()
projected = array([dot(V[:40], immatrix[i]-immean) for i in range(imnbr)])
cluster_num = 3
projected = whiten(projected)
centroids, distortion = kmeans(projected, cluster_num)
code, distance = vq(projected, centroids)
for k in range(cluster_num):
    ind = where(code==k)[0]
    figure()
    gray()
    for i in range(minimum(len(ind), 40)):
        subplot(4, 10, i+1)
        imshow(immatrix[ind[i]].reshape((25, 25)))
        axis('off')
show()

k=3でk平均法を用いたもの

#ブランチを作る関数
def divide_branch(data, branch, k):
    div = min(k, len(branch))
    if div<=1:
        return list(branch)
    centroids, distortion = kmeans(data[branch], k)
    code, distance = vq(data[branch], centroids)
    new_branch = []
    for i in range(k):
        ind = where(code==i)[0]
        if len(ind)==0:
            continue
        else:
            new_branch.append(divide_branch(data, branch[ind], k))
    return new_branch

#k=4で階層型k平均法を実行
tree = array([i for i in range(projected.shape[0])])
branches = divide_branch(projected, tree, 4)

#以下は描画のため
def get_depth(t):
    if len(t)<2:
        return 1
    else:
        return max([get_depth(tt) for tt in t])+1

def get_height(t):
    if (len(t)<2):
        return 1
    else:
        return sum([get_height(tt) for tt in t])

#ノードを描画する関数
def draw_node(node, draw, x, y, s, iml, im):
    if len(node)<1:
        return
    if len(node)==1:
        nodeim = Image.open(iml[node[0]])
        nodeim.thumbnail([20, 20])
        ns = nodeim.size
        im.paste(nodeim, [int(x), int(y-ns[1]//2), int(x+ns[0]), int(y+ns[1]-ns[1]//2)])
    else:
        ht = sum([get_height(n) for n in node])*20/2
        h1 = get_height(node[0])*20/2
        h2 = get_height(node[-1])*20/2
        top = y-ht
        bottom = y+ht
        draw.line((x, top+h1, x, bottom-h2), fill=(0, 0, 0))
        ll = 1*s
        y = top
        for i in range(len(node)):
            y += get_height(node[i])*20/2
            draw.line((x, y, x+ll, y), fill=(0, 0, 0))
            draw_node(node[i], draw, x+ll, y, s, imlist, im)
            y += get_height(node[i])*20/2

#ツリーを描画
def draw_dendrogram(node, iml, filename='kclusters.jpg'):
    rows = get_height(node)*20+40
    cols = 1200

    s = float(cols-150)/get_depth(node)

    im =  Image.new('RGB', (cols, rows), (255, 255, 255))
    draw = ImageDraw.Draw(im)

    draw.line((0, rows/2, 20, rows/2), fill=(0, 0, 0))
    draw_node(node, draw, 20, (rows/2), s, iml, im)
    im.save(filename)
    im.show()

#ツリーを出力
draw_dendrogram(branches, imlist, filename='k_fonts.jpg')

k=4 なので４つ以下の組にわけられています。あってそうですね。

6.4.2 Using the hierarchical k-means from the previous exercise, make a tree visualization (similar to the dendrogram for hierarchical clustering) that shows the average image for each cluster node. Tip: you can take the average PCA coefficients feature vector and use the PCA basis to synthesize an image for each feature vector.

演習問題１の結果を使って、各ノードの平均画像を表示するようにせよ

私の回答

問題のヒントとしてPCAのベクトルが使えるとありますので、それを利用します
https://github.com/moizumi99/CVBookExercise/blob/master/Chapter-6/CV%20Book%20Ch%206%20Exercise%202.ipynb

from PIL import Image, ImageDraw
from numpy import *
from pylab import *
import scipy.misc
from scipy.cluster.vq import *
import imtools
import pickle

imlist = imtools.get_imlist('selected_fontimages/')
imnbr = len(imlist)
with open('font_pca_modes.pkl', 'rb') as f:
    immean = pickle.load(f)
    V = pickle.load(f)

immatrix = array([array(Image.open(im)).flatten() for im in imlist], 'f')
immean = immean.flatten()
projected = array([dot(V[:40], immatrix[i]-immean) for i in range(imnbr)])
cluster_num = 3
projected = whiten(projected)
centroids, distortion = kmeans(projected, cluster_num)
code, distance = vq(projected, centroids)

def divide_branch_with_center(data, branch, k):
    div = min(k, len(branch))
    if div<=1:
        return list(branch)
    centroids, distortion = kmeans(data[branch], k)
    code, distance = vq(data[branch], centroids)
    new_branch = []
    for i in range(k):
        ind = where(code==i)[0]
        if len(ind)==0:
            continue
        else:
　　　　　　　#演習６．４．１とここが違う
            new_branch.append((centroids[i], distance[i], divide_branch_with_center(data, branch[ind], k)))
    return new_branch

#ツリーの作成
tree = array([i for i in range(projected.shape[0])])
branches = ([0 for i in range(40)], 0, divide_branch_with_center(projected, tree, 4))

#以下はツリーの描画用
def get_depth(t):
    if len(t[2])<2:
        return 1
    else:
        return max([get_depth(tt) for tt in t[2]])+1

def get_height(t):
    if (len(t[2])<2):
        return 1
    else:
        return sum([get_height(tt) for tt in t[2]])

def draw_average(center, x, y, im):
    c = center/np.linalg.norm(center)
    avim = dot((V[:40]).T, c)
    avim = 255*(avim-min(avim))/(max(avim)-min(avim)+1e-6)
    avim = avim.reshape(25, 25)
    avim[avim<0] = 0
    avim[avim>255] = 255
    avim = Image.fromarray(avim)
    avim.thumbnail([20, 20])
    ns = avim.size
    im.paste(avim, [int(x), int(y-ns[1]//2), int(x+ns[0]), int(y+ns[1]-ns[1]//2)])

def draw_node(node, draw, x, y, s, iml, im):
    if len(node[2])<1:
        return
    if len(node[2])==1:
        nodeim = Image.open(iml[node[2][0]])
        nodeim.thumbnail([20, 20])
        ns = nodeim.size
        im.paste(nodeim, [int(x), int(y-ns[1]//2), int(x+ns[0]), int(y+ns[1]-ns[1]//2)])
    else:
        ht = sum([get_height(n) for n in node[2]])*20/2
        h1 = get_height(node[2][0])*20/2
        h2 = get_height(node[2][-1])*20/2
        top = y-ht
        bottom = y+ht
        draw.line((x, top+h1, x, bottom-h2), fill=(0, 0, 0))
        y = top
        for i in range(len(node[2])):
            ll = node[2][i][1]/8*s
            y += get_height(node[2][i])*20/2
            xx = x + ll + s/4
            draw.line((x, y, xx, y), fill=(0, 0, 0))
            if len(node[2][i][2])>1:
                draw_average(node[2][i][0], xx, y, im)
                xx = xx+20
            draw.line((xx, y, xx+s/4, y), fill=(0, 0, 0))
            xx = xx+s/4
            draw_node(node[2][i], draw, xx, y, s, imlist, im)
            y += get_height(node[2][i])*20/2

def draw_dendrogram(node, iml, filename='kclusters.jpg'):
    rows = get_height(node)*20+40
    cols = 1200

    s = float(cols-150)/get_depth(node)

    im =  Image.new('RGB', (cols, rows), (255, 255, 255))
    draw = ImageDraw.Draw(im)

    x = 0
    y = rows/2
    avim = Image.fromarray(immean.reshape(25, 25))
    avim.thumbnail([20, 20])
    ns = avim.size
    im.paste(avim, [int(x), int(y-ns[1]//2), int(x+ns[0]), int(y+ns[1]-ns[1]//2)])
    draw.line((x+20, y, x+40, y), fill=(0, 0, 0))
    draw_node(node, draw, x+40, (rows/2), s, iml, im)
    im.save(filename)
    im.show()

draw_dendrogram(branches, imlist, filename='k_fonts.jpg')

あれ？なんだか変な枠がついてますね。正規化がうまくいっていないのでしょうか。それ以外はあっていそうです。

6.4.3 By modifying the class used for hierarchical clustering to include the number of images below the node you have a simple and fast way of finding similar (tight) groups of a given size. Implement this small change and try it out on some real data. How does it perform?

階層クラスタリングの際、ノードに含まれる画像の数を記録するようにし、任意の数の類似画像を高速化せよ

私の回答

準備としてまず、似た画像を集めます。今回はFlickrからsunsetをキーワードに画像をあつめ、flickr-sunsetsというフォルダにおいておきました。

次にhcluster.pyを変更
以下の行をclass ClusterNode(object)の定義に追加

hcluster.py

    def extract_clusters_by_num(self, num):
        """ Extract list of sub-tree clusters from
        hcluster tree with count<num. """

        if self.count < num:
            return [self]
        return self.left.extract_clusters_by_num(num) + self.right.extract_clusters_by_num(num)

また、def hcluster内の次の行を

hcluster.py

     new_node = ClusterNode(new_vec,left=ni,right=nj,
                          distance=closest)

以下の様に変更

hcluster.py

        new_node = ClusterNode(new_vec, left=ni, right=nj, distance=closest, count=cnt)

この上で実行

from PIL import Image
from numpy import *
from pylab import *
import scipy.misc
from scipy.cluster.vq import *
from scipy.misc import imresize
import os
import hcluster

path = 'flickr-sunsets/'
imlist = [os.path.join(path, f) for f in os.listdir(path) if f.endswith('.jpg')]
features = zeros([len(imlist), 512])
for i, f in enumerate(imlist):
    im = array(Image.open(f))

    h, edges = histogramdd(im.reshape(-1, 3), 8, normed=True, range=[(0,255), (0, 255), (0, 255)])

    features[i] = h.flatten()
tree = hcluster.hcluster(features)

#新しい関数を実行して８個以下のクラスターを取得
clusters = tree.extract_clusters_by_num(8)

for c in clusters:
    elements = c.get_cluster_elements()
    nbr_elements = len(elements)
    if nbr_elements>3:
        figure(figsize=(8, 8))
        for p in range(minimum(nbr_elements, 20)):
            subplot(4, 8, p+1)
            im = array(Image.open(imlist[elements[p]]))
            imshow(im)
            axis('off')
show()

6.4.4 Experiment with using single and complete linking for building the hierarchical cluster tree. How do the resulting clusters differ?

階層クラスタリングで単リンク法と完全リンク法を使え

私の回答

まず、hcluster.pyに以下の変更を加えます

hcluster.py

from itertools import product

#中略

def L2dist(v1, v2):
    av1 = average(v1, axis=0)
    av2 = average(v2, axis=0)
    return sqrt(sum((av1-av2)**2))

def L2dist_single(v1, v2):
    return min([sqrt(sum((a-b)**2)) for a, b in product(v1, v2)])

def L2dist_complete(v1, v2):
    return max([sqrt(sum((a-b)**2)) for a, b in product(v1, v2)])

def L1dist(v1, v2):
    av1 = average(v1)
    av2 = average(v2)
    return sum(abs(av1-av2))

def L1dist_single(v1, v2):
    return min([sum(abs(a-b)) for a, b in product(v1, v2)])

def L1dist_complete(v1, v2):
    return max([sum(abs(a-b)) for a, b in product(v1, v2)])

単リンク法、完全リンク法、それぞれの距離関数をつくります。

次に計算してみます。

from PIL import Image
from numpy import *
from pylab import *
import scipy.misc
from scipy.cluster.vq import *
from scipy.misc import imresize
import os
import hcluster

path = 'alamo/'
imlist = [os.path.join(path, f) for f in os.listdir(path) if f.endswith('.jpg')]

features = zeros([len(imlist), 512])
for i, f in enumerate(imlist):
    im = array(Image.open(f))

    h, edges = histogramdd(im.reshape(-1, 3), 8, normed=True, range=[(0,255), (0, 255), (0, 255)])

    features[i] = h.flatten()

#以下通常の平均リンク
tree = hcluster.hcluster(features, hcluster_single.L2dist)
clusters = tree.extract_clusters(0.24*tree.distance)

for c in clusters:
    elements = c.get_cluster_elements()
    nbr_elements = len(elements)
    if nbr_elements>2:
        figure(figsize=(8, 8))
        for p in range(minimum(nbr_elements, 20)):
            subplot(4, 5, p+1)
            im = array(Image.open(imlist[elements[p]]))
            imshow(im)
            axis('off')
show()
hcluster.draw_dendrogram(tree, imlist, filename='alamo_average.pdf')

#単リンク法
tree = hcluster.hcluster(features, hcluster_single.L2dist_single)
clusters = tree.extract_clusters(0.24*tree.distance)

for c in clusters:
    elements = c.get_cluster_elements()
    nbr_elements = len(elements)
    if nbr_elements>2:
        figure(figsize=(8, 8))
        for p in range(minimum(nbr_elements, 20)):
            subplot(4, 5, p+1)
            im = array(Image.open(imlist[elements[p]]))
            imshow(im)
            axis('off')
show()
hcluster.draw_dendrogram(tree, imlist, filename='alamo_single.pdf')

#完全リンク法
tree = hcluster.hcluster(features, hcluster_single.L2dist_complete)
clusters = tree.extract_clusters(0.18*tree.distance)

for c in clusters:
    elements = c.get_cluster_elements()
    nbr_elements = len(elements)
    if nbr_elements>2:
        figure(figsize=(8, 8))
        for p in range(minimum(nbr_elements, 20)):
            subplot(4, 5, p+1)
            im = array(Image.open(imlist[elements[p]]))
            imshow(im)
            axis('off')
show()
hcluster.draw_dendrogram(tree, imlist, filename='alamo_complete.pdf')

クラスター化される画像の傾向が多少違うようですが、どう解釈してよいのかよくわかりません

6.4.5 In some spectral clustering algorithms the matrix D1S

is used instead of L. Try replacing the Laplacian matrix with this and apply this on a few different data sets.
D^-1Sを使ってスペクトラルクラスタリングを実装せよ

私の回答

'golden gate bridge'というキーワードでflickrで検索した画像を使いました。

from PIL import Image
from numpy import *
from pylab import *
import scipy.misc
from scipy.cluster.vq import *
import imtools
import pickle
import sift
import os

imlist = imtools.get_imlist('./goldengatebridge/')
nbr_images = len(imlist)
matchscores = zeros((nbr_images, nbr_images))
featlist = []
for imname in imlist:
    name, ext = os.path.splitext(imname)
    featname = name+'.sift'
    sift.process_image(imname, featname)
    featlist.append(featname)

cnt = 0
for i in range(nbr_images):
    for j in range(i, nbr_images):
        cnt += 1
        print cnt, "/ ", (nbr_images+1)*nbr_images/2
        print 'comparing', imlist[i], imlist[j]

        try:
            l1, d1 = sift.read_features_from_file(featlist[i])
            l2, d2 = sift.read_features_from_file(featlist[j])
            matches = sift.match_twosided(d1, d2)
        except(IndexError):
            matches = []

        nbr_matches = sum(matches > 0)
        print 'number of matches = ', nbr_matches
        matchscores[i, j] = nbr_matches

for i in range(nbr_images):
    for j in range(i+1, nbr_images):
        matchscores[j, i] = matchscores[i, j]

print matchscores

[[ 462. 0. 0. ..., 0. 0. 0.]
[ 0. 176. 0. ..., 1. 0. 0.]
[ 0. 0. 159. ..., 0. 0. 0.]
...,
[ 0. 1. 0. ..., 410. 0. 1.]
[ 0. 0. 0. ..., 0. 495. 0.]
[ 0. 0. 0. ..., 1. 0. 283.]]

#まずは通常のラプラシアンで
n = len(imlist)
S = 1/(matchscores+1e-6)

rowsum = sum(S, axis=0)
D = diag(1/sqrt(rowsum))
I = identity(n)
L = I - dot(D, dot(S, D))
U, sigma, V = linalg.svd(L)

k = 4
features = array(V[:k]).T
features = whiten(features)
centroids, distortion = kmeans(features, k)
code, distance = vq(features, centroids)

for c in range(k):
    ind = where(code==c)[0]
    figure(figsize=(16, 8))
    for i in range(minimum(len(ind), 39)):
        im = Image.open(imlist[ind[i]])
        subplot(4, 10, i+1)
        gray()
        imshow(array(im))
        axis('equal')
        axis('off')
show()

ラプラシアンによるスペクトラルクラスタリング

#次にD^-1Sによるもの
n = len(imlist)
S = 1/(matchscores+1e-6)
rowsum = sum(S, axis=0)
D = diag(1/rowsum)
I = identity(n)
L = dot(D, S)
U, sigma, V = linalg.svd(L)
k = 4
features = array(V[:k]).T
features = whiten(features)
centroids, distortion = kmeans(features, k)
code, distance = vq(features, centroids)

for c in range(k):
    ind = where(code==c)[0]
    figure(figsize=(16, 8))
    for i in range(minimum(len(ind), 39)):
        im = Image.open(imlist[ind[i]])
        subplot(4, 10, i+1)
        gray()
        imshow(array(im))
        axis('equal')
        axis('off')
show()

D^-1Sによるスペクトラルクラスタリング

ラプラシアンの結果のほうが良いような気がしますね。

6.4.6 Download some image collections from Flickr searches with different tags. Extract the RGB histogram like you did for the sunset images. Cluster the images using one of the methods from this chapter. Can you separate the classes with the clusters?

フリッカーからいくつかイメージのセットを違うタグであつめて、RGBヒストグラムを計算し、クラスター化せよ

私の回答

とりあえず、またゴールデン・ゲート・ブリッジ画像で

from PIL import Image
from numpy import *
from pylab import *
import scipy.misc
from scipy.cluster.vq import *
from scipy.misc import imresize
import os
import hcluster

path = 'goldengatebridge/'
imlist = [os.path.join(path, f) for f in os.listdir(path) if f.endswith('.jpg')]

#sunsetと同じように
features = zeros([len(imlist), 512])
for i, f in enumerate(imlist):
    im = array(Image.open(f))

    h, edges = histogramdd(im.reshape(-1, 3), 8, normed=True, range=[(0,255), (0, 255), (0, 255)])

    features[i] = h.flatten()

tree = hcluster.hcluster(features)
clusters = tree.extract_clusters(0.6*tree.distance)

for c in clusters:
    elements = c.get_cluster_elements()
    nbr_elements = len(elements)
    if nbr_elements>3:
        figure(figsize=(8, 8))
        for p in range(minimum(nbr_elements, 20)):
            subplot(4, 5, p+1)
            im = array(Image.open(imlist[elements[p]]))
            imshow(im)
            axis('off')
show()

#スペクトラルクラスタリングを適用してみる
n = len(features)
S = array([[sqrt(sum((features[i]-features[j])**2)) for i in range(n)] for j in range(n)], 'f')
rowsum = sum(S, axis=0)
D = diag(1/sqrt(rowsum))
I = identity(n)
L = I - dot(D, dot(S, D))
U, sigma, V = linalg.svd(L)

k = 5
f = array(V[:k]).T
f = whiten(f)
centroids, distortion = kmeans(f, k)
code, distance = vq(f, centroids)

#結果の描画
for c in range(k):
    ind = where(code==c)[0]
    figure()
    for i in range(minimum(len(ind), 39)):
        im = Image.open(imlist[ind[i]])
        subplot(4, 10, i+1)
        gray()
        imshow(array(im))
        axis('equal')
        axis('off')
show()