图像聚类

1. 6.1 K-Means聚类
   1. 6.1.1 SciPy聚类包
   2. 6.1.2 图像聚类
   3. 6.1.3 在主成分上可视化图像
   4. 6.1.4 像素聚类
2. 6.2 层次聚类
   1. 6.2.1 图像聚类
3. 6.3 谱聚类

这一章会介绍几种聚类方法，并就怎么使用它们对图像进行聚类找出相似的图像组进行说明。聚类可以用于识别，划分图像数据集、组织导航等。同时，我们也会用聚类相似的图像进行可视化。

6.1 K-Means聚类

K-means是一种非常简单的聚类算法，它能够将输入数据划分成k个簇。关于K-means聚类算法的介绍可以参阅中译本。

6.1.1 SciPy聚类包

尽管K-means聚类算法很容易实现，但我们没必要自己去实现。SciPy矢量量化包sci.cluter.vq中有k-means的实现。这里我们演示怎样使用它。

我们以2维示例样本数据进行说明：

1. # coding=utf-8

2. """

3. Function: figure 6.1

4. An example of k-means clustering of 2D points

5. """

6. from pylab import *

7. from scipy.cluster.vq import *

8.  

9. # 添加中文字体支持

10. from matplotlib.font_manager import FontProperties

11. font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)

12.  

13. class1 = 1.5 * randn(100, 2)

14. class2 = randn(100, 2) + array([5, 5])

15. features = vstack((class1, class2))

16. centroids, variance = kmeans(features, 2)

17. code, distance = vq(features, centroids)

18. figure()

19. ndx = where(code == 0)[0]

20. plot(features[ndx, 0], features[ndx, 1], '*')

21. ndx = where(code == 1)[0]

22. plot(features[ndx, 0], features[ndx, 1], 'r.')

23. plot(centroids[:, 0], centroids[:, 1], 'go')

24.  

25. title(u'2维数据点聚类', fontproperties=font)

26. axis('off')

27. show()

上面代码中where()函数给出每类的索引。运行上面代码，可得到原书P129页图6-1，即：

6.1.2 图像聚类

现在我们用k-means对原书14页的图像进行聚类，文件selectedfontimages.zip包含了66张字体图像。对于每一张图像，我们用在前40个主成分上投影后的系数作为特征向量。下面为对其进行聚类的代码：

1. # -*- coding: utf-8 -*-

2. from PCV.tools import imtools

3. import pickle

4. from scipy import *

5. from pylab import *

6. from PIL import Image

7. from scipy.cluster.vq import *

8. from PCV.tools import pca

9.  

10. # Uses sparse pca codepath.

11. imlist = imtools.get_imlist('../data/selectedfontimages/a_selected_thumbs/')

12.  

13. # 获取图像列表和他们的尺寸

14. im = array(Image.open(imlist[0])) # open one image to get the size

15. m, n = im.shape[:2] # get the size of the images

16. imnbr = len(imlist) # get the number of images

17. print "The number of images is %d" % imnbr

18.  

19. # Create matrix to store all flattened images

20. immatrix = array([array(Image.open(imname)).flatten() for imname in imlist], 'f')

21.  

22. # PCA降维

23. V, S, immean = pca.pca(immatrix)

24.  

25. # 保存均值和主成分

26. #f = open('./a_pca_modes.pkl', 'wb')

27. f = open('./a_pca_modes.pkl', 'wb')

28. pickle.dump(immean,f)

29. pickle.dump(V,f)

30. f.close()

31.  
32.  

33. # get list of images

34. imlist = imtools.get_imlist('../data/selectedfontimages/a_selected_thumbs/')

35. imnbr = len(imlist)

36.  

37. # load model file

38. with open('../data/selectedfontimages/a_pca_modes.pkl','rb') as f:

39. immean = pickle.load(f)

40. V = pickle.load(f)

41. # create matrix to store all flattened images

42. immatrix = array([array(Image.open(im)).flatten() for im in imlist],'f')

43.  

44. # project on the 40 first PCs

45. immean = immean.flatten()

46. projected = array([dot(V[:40],immatrix[i]-immean) for i in range(imnbr)])

47.  

48. # k-means

49. projected = whiten(projected)

50. centroids,distortion = kmeans(projected,4)

51. code,distance = vq(projected,centroids)

52.  

53. # plot clusters

54. for k in range(4):

55. ind = where(code==k)[0]

56. figure()

57. gray()

58. for i in range(minimum(len(ind),40)):

59. subplot(4,10,i+1)

60. imshow(immatrix[ind[i]].reshape((25,25)))

61. axis('off')

62. show()

运行上面代码，可得到下面的聚类结果：注：这里的结果译者截的是原书上的结果，上面代码实际运行出来的结果可能跟上面有出入。

6.1.3 在主成分上可视化图像

1. # -*- coding: utf-8 -*-

2. from PCV.tools import imtools, pca

3. from PIL import Image, ImageDraw

4. from pylab import *

5. from PCV.clustering import hcluster

6.  

7. imlist = imtools.get_imlist('../data/selectedfontimages/a_selected_thumbs')

8. imnbr = len(imlist)

9.  

10. # Load images, run PCA.

11. immatrix = array([array(Image.open(im)).flatten() for im in imlist], 'f')

12. V, S, immean = pca.pca(immatrix)

13.  

14. # Project on 2 PCs.

15. projected = array([dot(V[[0, 1]], immatrix[i] - immean) for i in range(imnbr)]) # P131 Fig6-3左图

16. #projected = array([dot(V[[1, 2]], immatrix[i] - immean) for i in range(imnbr)]) # P131 Fig6-3右图

17.  

18. # height and width

19. h, w = 1200, 1200

20.  

21. # create a new image with a white background

22. img = Image.new('RGB', (w, h), (255, 255, 255))

23. draw = ImageDraw.Draw(img)

24.  

25. # draw axis

26. draw.line((0, h/2, w, h/2), fill=(255, 0, 0))

27. draw.line((w/2, 0, w/2, h), fill=(255, 0, 0))

28.  

29. # scale coordinates to fit

30. scale = abs(projected).max(0)

31. scaled = floor(array([(p/scale) * (w/2 - 20, h/2 - 20) + (w/2, h/2)

32. for p in projected])).astype(int)

33.  

34. # paste thumbnail of each image

35. for i in range(imnbr):

36. nodeim = Image.open(imlist[i])

37. nodeim.thumbnail((25, 25))

38. ns = nodeim.size

39. box = (scaled[i][0] - ns[0] // 2, scaled[i][1] - ns[1] // 2,

40. scaled[i][0] + ns[0] // 2 + 1, scaled[i][1] + ns[1] // 2 + 1)

41. img.paste(nodeim, box)

42.  

43. tree = hcluster.hcluster(projected)

44. hcluster.draw_dendrogram(tree,imlist,filename='fonts.png')

45.  

46. figure()

47. imshow(img)

48. axis('off')

49. img.save('../images/ch06/pca_font.png')

50. show()

运行上面代码，可画出原书P131图6-3中的实例结果。

6.1.4 像素聚类

在结束这节前，我们看一个对像素进行聚类而不是对所有的图像进行聚类的例子。将图像区域归并成“有意义的”组件称为图像分割。在第九章会将其单独列为一个主题。在像素级水平进行聚类除了可以用在一些很简单的图像，在其他图像上进行聚类是没有意义的。这里，我们将k-means应用到RGB颜色值上，关于分割问题会在第九章第二节会给出分割的方法。下面是对两幅图像进行像素聚类的例子(注：译者对原书中的代码做了调整)：

1. # -*- coding: utf-8 -*-

2. """

3. Function: figure 6.4

4. Clustering of pixels based on their color value using k-means.

5. """

6. from scipy.cluster.vq import *

7. from scipy.misc import imresize

8. from pylab import *

9. import Image

10.  

11. # 添加中文字体支持

12. from matplotlib.font_manager import FontProperties

13. font = FontProperties(fname=r"c:\windows\fonts\SimSun.ttc", size=14)

14.  

15. def clusterpixels(infile, k, steps):

16. im = array(Image.open(infile))

17. dx = im.shape[0] / steps

18. dy = im.shape[1] / steps

19. # compute color features for each region

20. features = []

21. for x in range(steps):

22. for y in range(steps):

23. R = mean(im[x * dx:(x + 1) * dx, y * dy:(y + 1) * dy, 0])

24. G = mean(im[x * dx:(x + 1) * dx, y * dy:(y + 1) * dy, 1])

25. B = mean(im[x * dx:(x + 1) * dx, y * dy:(y + 1) * dy, 2])

26. features.append([R, G, B])

27. features = array(features, 'f') # make into array

28. # 聚类， k是聚类数目

29. centroids, variance = kmeans(features, k)

30. code, distance = vq(features, centroids)

31. # create image with cluster labels

32. codeim = code.reshape(steps, steps)

33. codeim = imresize(codeim, im.shape[:2], 'nearest')

34. return codeim

35.  

36. k=3

37. infile_empire = '../data/empire.jpg'

38. im_empire = array(Image.open(infile_empire))

39. infile_boy_on_hill = '../data/boy_on_hill.jpg'

40. im_boy_on_hill = array(Image.open(infile_boy_on_hill))

41. steps = (50, 100) # image is divided in steps*steps region

42. print steps[0], steps[-1]

43.  

44. #显示原图empire.jpg

45. figure()

46. subplot(231)

47. title(u'原图', fontproperties=font)

48. axis('off')

49. imshow(im_empire)

50.  

51. # 用50*50的块对empire.jpg的像素进行聚类

52. codeim= clusterpixels(infile_empire, k, steps[0])

53. subplot(232)

54. title(u'k=3,steps=50', fontproperties=font)

55. #ax1.set_title('Image')

56. axis('off')

57. imshow(codeim)

58.  

59. # 用100*100的块对empire.jpg的像素进行聚类

60. codeim= clusterpixels(infile_empire, k, steps[-1])

61. ax1 = subplot(233)

62. title(u'k=3,steps=100', fontproperties=font)

63. #ax1.set_title('Image')

64. axis('off')

65. imshow(codeim)

66.  

67. #显示原图empire.jpg

68. subplot(234)

69. title(u'原图', fontproperties=font)

70. axis('off')

71. imshow(im_boy_on_hill)

72.  

73. # 用50*50的块对empire.jpg的像素进行聚类

74. codeim= clusterpixels(infile_boy_on_hill, k, steps[0])

75. subplot(235)

76. title(u'k=3,steps=50', fontproperties=font)

77. #ax1.set_title('Image')

78. axis('off')

79. imshow(codeim)

80.  

81. # 用100*100的块对empire.jpg的像素进行聚类

82. codeim= clusterpixels(infile_boy_on_hill, k, steps[-1])

83. subplot(236)

84. title(u'k=3，steps=100', fontproperties=font)

85. axis('off')

86. imshow(codeim)

87.  

88. show()

上面代码中，先载入一幅图像，然后用一个steps*steps的方块在原图中滑动，对窗口中的图像值求和取平均，将它下采样到一个较低的分辨率，然后对这些区域用k-means进行聚类。运行上面代码，即可得出原书P133页图6-4中的图。

6.2 层次聚类

层次聚类(或称凝聚聚类)是另一种简单但有效的聚类算法。下面我们我们通过一个简单的实例看看层次聚类是怎样进行的。

1. from pylab import *

2. from PCV.clustering import hcluster

3.  

4. class1 = 1.5 * randn(100,2)

5. class2 = randn(100,2) + array([5,5])

6. features = vstack((class1,class2))

7.  

8. tree = hcluster.hcluster(features)

9. clusters = tree.extract_clusters(5)

10. print 'number of clusters', len(clusters)

11. for c in clusters:

12. print c.get_cluster_elements()

上面代码首先创建一些2维数据点，然后对这些数据点聚类，用一些阈值提取列表中的聚类后的簇群，并将它们打印出来，译者在自己的笔记本上打印出的结果为：

number of clusters 2
[197, 107, 176, 123, 173, 189, 154, 136, 183, 113, 109, 199, 178, 129, 163, 100, 148, 111, 143, 118, 162, 169, 138, 182, 193, 116, 134, 198, 184, 181, 131, 166, 127, 185, 161, 171, 152, 157, 112, 186, 128, 156, 108, 158, 120, 174, 102, 137, 117, 194, 159, 105, 155, 132, 188, 125, 180, 151, 192, 164, 195, 126, 103, 196, 179, 146, 147, 135, 139, 110, 140, 106, 104, 115, 149, 190, 170, 172, 121, 145, 114, 150, 119, 142, 122, 144, 160, 187, 153, 167, 130, 133, 165, 191, 175, 177, 101, 141, 124, 168]
[0, 39, 32, 87, 40, 48, 28, 8, 26, 12, 94, 5, 1, 61, 24, 59, 83, 10, 99, 50, 23, 58, 51, 16, 71, 25, 11, 37, 22, 46, 60, 86, 65, 2, 21, 4, 41, 72, 80, 84, 33, 56, 75, 77, 29, 85, 93, 7, 73, 6, 82, 36, 49, 98, 79, 43, 91, 14, 47, 63, 3, 97, 35, 18, 44, 30, 13, 67, 62, 20, 57, 89, 88, 9, 54, 19, 15, 92, 38, 64, 45, 70, 52, 95, 69, 96, 42, 53, 27, 66, 90, 81, 31, 34, 74, 76, 17, 78, 55, 68]

6.2.1 图像聚类

1. # -*- coding: utf-8 -*-

2. import os

3. import Image

4. from PCV.clustering import hcluster

5. from matplotlib.pyplot import *

6. from numpy import *

7.  

8. # create a list of images

9. path = '../data/sunsets/flickr-sunsets-small/'

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

11. # extract feature vector (8 bins per color channel)

12. features = zeros([len(imlist), 512])

13. for i, f in enumerate(imlist):

14. im = array(Image.open(f))

15. # multi-dimensional histogram

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

17. features[i] = h.flatten()

18. tree = hcluster.hcluster(features)

19.  

20. # visualize clusters with some (arbitrary) threshold

21. clusters = tree.extract_clusters(0.23 * tree.distance)

22. # plot images for clusters with more than 3 elements

23. for c in clusters:

24. elements = c.get_cluster_elements()

25. nbr_elements = len(elements)

26. if nbr_elements > 3:

27. figure()

28. for p in range(minimum(nbr_elements,20)):

29. subplot(4, 5, p + 1)

30. im = array(Image.open(imlist[elements[p]]))

31. imshow(im)

32. axis('off')

33. show()

34.  

35. hcluster.draw_dendrogram(tree,imlist,filename='sunset.pdf')

运行上面代码，可得原书P140图6-6。同时会在上面脚本文件所在的文件夹下生成层次聚类后的簇群树：我们对前面字体图像同样创建一个树，正如前面在主成分可视化图像中，我们添加了下面代码：

tree = hcluster.hcluster(projected)
hcluster.draw_dendrogram(tree,imlist,filename='fonts.png')

运行添加上面两行代码后前面的例子，可得对字体进行层次聚类后的簇群树：

6.3 谱聚类

谱聚类是另一种不同于k-means和层次聚类的聚类算法。关于谱聚类的原理，可以参阅中译本。这里，我们用原来k-means实例中用到的字体图像。

1. # -*- coding: utf-8 -*-

2. from PCV.tools import imtools, pca

3. from PIL import Image, ImageDraw

4. from pylab import *

5. from scipy.cluster.vq import *

6.  

7. imlist = imtools.get_imlist('../data/selectedfontimages/a_selected_thumbs')

8. imnbr = len(imlist)

9.  

10. # Load images, run PCA.

11. immatrix = array([array(Image.open(im)).flatten() for im in imlist], 'f')

12. V, S, immean = pca.pca(immatrix)

13.  

14. # Project on 2 PCs.

15. projected = array([dot(V[[0, 1]], immatrix[i] - immean) for i in range(imnbr)]) # P131 Fig6-3左图

16. #projected = array([dot(V[[1, 2]], immatrix[i] - immean) for i in range(imnbr)]) # P131 Fig6-3右图

17.  

18. n = len(projected)

19. # compute distance matrix

20. S = array([[ sqrt(sum((projected[i]-projected[j])**2))

21. for i in range(n) ] for j in range(n)], 'f')

22. # create Laplacian matrix

23. rowsum = sum(S,axis=0)

24. D = diag(1 / sqrt(rowsum))

25. I = identity(n)

26. L = I - dot(D,dot(S,D))

27. # compute eigenvectors of L

28. U,sigma,V = linalg.svd(L)

29. k = 5

30. # create feature vector from k first eigenvectors

31. # by stacking eigenvectors as columns

32. features = array(V[:k]).T

33. # k-means

34. features = whiten(features)

35. centroids,distortion = kmeans(features,k)

36. code,distance = vq(features,centroids)

37. # plot clusters

38. for c in range(k):

39. ind = where(code==c)[0]

40. figure()

41. gray()

42. for i in range(minimum(len(ind),39)):

43. im = Image.open(imlist[ind[i]])

44. subplot(4,10,i+1)

45. imshow(array(im))

46. axis('equal')

47. axis('off')

48. show()

上面我们在前个特征向量上计算标准的k-means。下面是运行上面代码的结果：注意，由于在k-means阶段会给出不同的聚类结果，所以你运行上面代码出来的结果可能跟译者的是不一样的。

同样，我们可以在不知道特征向量或是没有严格相似性定义的情况下进行谱聚类。原书44页的位置地理图像是通过它们之间有多少局部描述子匹配相连接的。48页的相似性矩阵中的元素是为规范化的匹配特征点数。我们同样可以对其进行谱聚类，完整的代码如下：

1. # -*- coding: utf-8 -*-

2. from PCV.tools import imtools, pca

3. from PIL import Image, ImageDraw

4. from PCV.localdescriptors import sift

5. from pylab import *

6. import glob

7. from scipy.cluster.vq import *

8.  
9.  

10. #download_path = "panoimages" # set this to the path where you downloaded the panoramio images

11. #path = "/FULLPATH/panoimages/" # path to save thumbnails (pydot needs the full system path)

12.  

13. download_path = "F:/dropbox/Dropbox/translation/pcv-notebook/data/panoimages" # set this to the path where you downloaded the panoramio images

14. path = "F:/dropbox/Dropbox/translation/pcv-notebook/data/panoimages/" # path to save thumbnails (pydot needs the full system path)

15.  

16. # list of downloaded filenames

17. imlist = imtools.get_imlist('../data/panoimages/')

18. nbr_images = len(imlist)

19.  

20. # extract features

21. #featlist = [imname[:-3] + 'sift' for imname in imlist]

22. #for i, imname in enumerate(imlist):

23. # sift.process_image(imname, featlist[i])

24.  

25. featlist = glob.glob('../data/panoimages/*.sift')

26.  

27. matchscores = zeros((nbr_images, nbr_images))

28.  

29. for i in range(nbr_images):

30. for j in range(i, nbr_images): # only compute upper triangle

31. print 'comparing ', imlist[i], imlist[j]

32. l1, d1 = sift.read_features_from_file(featlist[i])

33. l2, d2 = sift.read_features_from_file(featlist[j])

34. matches = sift.match_twosided(d1, d2)

35. nbr_matches = sum(matches > 0)

36. print 'number of matches = ', nbr_matches

37. matchscores[i, j] = nbr_matches

38. print "The match scores is: \n", matchscores

39.  

40. # copy values

41. for i in range(nbr_images):

42. for j in range(i + 1, nbr_images): # no need to copy diagonal

43. matchscores[j, i] = matchscores[i, j]

44.  

45. n = len(imlist)

46. # load the similarity matrix and reformat

47. S = matchscores

48. S = 1 / (S + 1e-6)

49. # create Laplacian matrix

50. rowsum = sum(S,axis=0)

51. D = diag(1 / sqrt(rowsum))

52. I = identity(n)

53. L = I - dot(D,dot(S,D))

54. # compute eigenvectors of L

55. U,sigma,V = linalg.svd(L)

56. k = 2

57. # create feature vector from k first eigenvectors

58. # by stacking eigenvectors as columns

59. features = array(V[:k]).T

60. # k-means

61. features = whiten(features)

62. centroids,distortion = kmeans(features,k)

63. code,distance = vq(features,centroids)

64. # plot clusters

65. for c in range(k):

66. ind = where(code==c)[0]

67. figure()

68. gray()

69. for i in range(minimum(len(ind),39)):

70. im = Image.open(imlist[ind[i]])

71. subplot(5,4,i+1)

72. imshow(array(im))

73. axis('equal')

74. axis('off')

75. show()

改变聚类数目k，可以得到不同的结果。译者分别测试了原书中k=2和k=10的情况，运行结果如下： k=2k=10注：对于聚类后，图像小于或等于1的类，在上面没有显示。

from: http://yongyuan.name/pcvwithpython/chapter6.html