[TOC]

基本的图像操作和处理

PIL库简介,文档

PIL是python的图像处理库。下面的例子是一些基本用法。

  1. 打开图像。通过open得到一个图像对象

    from PIL import Image
    pil_im = Image.open('test.jpg')
    
  2. 图像转换用convert,转换成灰度图用L

    pil_im = Image.open('../data/empire.jpg').convert('L')
    
  3. 创建缩略图

    pil_im.thumbnail((200,200))
    
  4. 复制和粘贴图像区域,crop和paste

  5. 调整图片大小,resize

    out = pil_im.resize((128, 128))
    
  6. 旋转图像,rotate

    out = pil_im.rotate(45)
    
  7. 显示图片

    pil_im.show()
    

Matplotlib库,文档

使用其中的Pylab接口绘制图表、点、线、曲线等类库。

  1. 使用plot来绘制点或线。有三个参数,分别为x轴的点列表,y轴的点列表,样式。

    from pylab import *
    
    #打开图像并转为数组
    im = array(Image.open('../data/empire.jpg'))
    
    #在pylab中需要调用imshow绘制图像,参数为数组
    imshow(im)
    
    x = [100, 100, 400, 400]
    y = [200, 500, 200, 500]
    
    #三个参数,绘制红色的x点
    plot(x, y, 'rx')
    
    #这两个参数在绘制线,分别对应四个点中的三条线
    plot(x, y)
    show()
    

    图像如下:

  2. 上图有坐标,可以使用axis(‘off’)来关闭

  3. 图像轮廓contour和直方图hist

    #得到一个灰度图的数组,是个二维数组[[1,2,3...],[3,2,4,...],...]
    im = array(Image.open('../data/empire.jpg').convert('L'))
    
    #新建图片,同一个脚本可以调用多次。展示时,关闭一个再展示下一个。
    figure()
    
    #不使用颜色信息
    gray()
    
    #显示轮廓信息
    contour(im, origin='image')
    
    #将刻度和比例设置为相同
    axis('equal')
    axis('off')
    
    figure()
    
    #flatten将二维数组转为一维数组
    imf = im.flatten()
    #print(im)
    #print(imf)
    
    #绘制直方图,只接受一维数组
    hist(imf, 128)
    show()
    

    得到图像:

NumPy库,文档

是一个科学计算工具包。在py中一般使用列表来表示数组,但是列表可以表示任何对象,所以[1,2,3]中其实保存的是对象的指针而不是数值,比较耗资源。所以在NumPy中提供了两种基本对象,ndarray用来存储单一数据类型的多维数组,ufunc为处理数组的函数。

  1. 图像的数组表示。使用NumPy的array方法保存。

    正常图像转数组,如下:

    [[[ 88 133 188]
     [ 86 131 186]
     [ 86 131 186]
     ...,
     [247 248 250]
     [246 247 249]
     [246 247 249]]
    
     [[ 89 134 189]
     [ 88 133 188]
     [ 88 133 188]
     ...,
     [246 247 249]
     [246 247 249]
     [246 247 249]]
    
     [[ 88 135 189]
     [ 87 134 188]
     [ 86 133 187]
     ...,
     [245 246 250]
     [245 246 250]
     [246 247 251]]
    
     ...,
     [[130 177 231]
     [131 178 232]
     [131 178 232]
     ...,
     [ 35 40 36]
     [ 36 41 37]
     [ 35 41 37]]
    
     [[130 177 231]
     [131 178 232]
     [131 178 232]
     ...,
     [ 36 41 37]
     [ 39 44 40]
     [ 37 43 39]]
    
     [[128 178 231]
     [130 180 233]
     [128 178 231]
     ...,
     [ 37 42 38]
     [ 36 42 38]
     [ 35 41 37]]]
    

    是个三维数组,包含行、列,以及RGB三个颜色通道。
    如果是灰度图像,如下:

    [[125 123 123 ..., 247 246 246]
     [126 125 125 ..., 246 246 246]
     [127 126 125 ..., 246 246 247]
     ...,
     [169 170 170 ..., 38 39 38]
     [169 170 170 ..., 39 42 40]
     [169 171 169 ..., 40 39 38]]
    

    是个二维数组。

  2. 灰度变换

    im = array(Image.open('../data/empire.jpg').convert('L'))
    
    #反相处理
    im2 = 255 - im
    im3 = (100.0/255) * im + 100
    
    #求平方后得到的图像
    im4 = 255.0 * (im/255.0)**2
    imshow(im2)
    show()
    

    如上,得到im2的图像:

  3. 图像缩放。NumPy都是数组,实现缩放需要先转为图像对象缩放再转回来。

    def imresize(im, sz):
     pil_im = Image.fromarray(uint8(im))
     return array(pil_im.resize(sz))
    
  4. 直方图均衡化。是指将一幅图像的灰度直方图变平,使变换后的图像中每个灰度值的分布概率都相同。

    def histeq(im,nbr_bins=256):
     """ Histogram equalization of a grayscale image. """
    
    # get image histogram
    
     imhist,bins = histogram(im.flatten(),nbr_bins,normed=True)
     cdf = imhist.cumsum() # cumulative distribution function
     cdf = 255 * cdf / cdf[-1] # normalize
    
    # use linear interpolation of cdf to find new pixel values
    
     im2 = interp(im.flatten(),bins[:-1],cdf)
     return im2.reshape(im.shape), cdf
    

    处理过的图像可以见对比:

    均衡化后,对比度增强了,细节也清晰了。

  5. 图像平均。是减少图像噪音的简单方式。一般是简单相加然后处理图像数目来计算平均图像。把多张相同尺寸图片合成一张。

  6. 图像的主成分分析(PCA)。

    比如一个100*100的灰度图像也有10000维,一兆像素则有百万维,所以我们需要降维操作。PCA是一个非常有用的降维技巧,可以在尽可能少维数的前提下尽量多地保持训练数据的信息。

    进行PCA变换图像要转为一维向量,使用flatten()方法来压平,一个25*25的图像的灰度数组压平后,为一个625的一维数组。类似[[1,2,3],[4,5,6],[7,8,9]]压平为[1,2,3,4,5,6,7,8,9]。

    比如例子里我们要对2359个不同a字母图像进行PCA变换,我们对每个字母图片压平成一个一维数组,然后2359个图片组合起来成一个(2359,625)的大数组。然后对这个大数组进行PCA变换,返回投影矩阵,方差,均值。最后我们使用reshape函数重新转换为二维图像。最后能得到这2359张图片的平均图片。

    PCA算法:

    def pca(X):
     """ Principal Component Analysis
     input: X, matrix with training data stored as flattened arrays in rows
     return: projection matrix (with important dimensions first), variance and mean.
     """
    
    # get dimensions
    
     num_data,dim = X.shape
    
    # center data
    
     mean_X = X.mean(axis=0)
     X = X - mean_X
    
     if dim>num_data:
    
    # PCA - compact trick used
    
     M = dot(X,X.T) # covariance matrix
     e,EV = linalg.eigh(M) # eigenvalues and eigenvectors
     tmp = dot(X.T,EV).T # this is the compact trick
     V = tmp[::-1] # reverse since last eigenvectors are the ones we want
     S = sqrt(e)[::-1] # reverse since eigenvalues are in increasing order
     for i in range(V.shape[1]):
     V[:,i] /= S
     else:
    
    # PCA - SVD used
    
     U,S,V = linalg.svd(X)
     V = V[:num_data] # only makes sense to return the first num_data
    
    # return the projection matrix, the variance and the mean
    
     return V,S,mean_X
    

    返回投影矩阵(二维,625*625)、方差(一维,625)和均值(一维,625)。
    运行代码:

    imlist = imtools.get_imlist('../a_thumbs')
    im = array(Image.open(imlist[0])) # open one image to get the size
    m, n = im.shape[:2] # get the size of the images
    imnbr = len(imlist) # get the number of images
    print "The number of images is %d" % imnbr
    
    # Create matrix to store all flattened images
    
    immatrix = array([array(Image.open(imname)).flatten() for imname in imlist], 'f')
    
    V, S, immean = pca.pca(immatrix)
    
    # Show the images (mean and 7 first modes)
    
    # This gives figure 1-8 (p15) in the book.
    
    figure()
    gray()
    
    #在一张图上生成多个图像
    subplot(2, 4, 1)
    axis('off')
    imshow(immean.reshape(m, n))
    for i in range(7):
     subplot(2, 4, i+2)
     imshow(V[i].reshape(m, n))
     axis('off')
    show()
    

    得到图片:
    其中第一张为平均图像,其余为前7个模式(具有最大方差的方向模式)

  7. 使用pickle模块把对象转换成字符串,以便传输和保存。一般保存到.pkl文件中。

SciPy库,文档

SciPy是建立在NumPy基础上用于数值运算的开源工具包。

  1. 高斯模糊。一般在图像插值操作、兴趣点计算等操作用会用到高斯模糊。它使用scipy.ndimage.filters中的gaussian_filter方法。

  2. 图像导数。

    在很多应用中图像强度的变化是非常重要的信息。强度的变化可以用灰度图像_I_(彩色图像一般对每个颜色通道分别计算导数)的x和y方向导数_Iχ_和_Iγ_进行描述。

    使用Sobel导数滤波器计算导数图像

    im = array(Image.open('../data/empire.jpg').convert('L'))
    gray()
    
    subplot(1, 4, 1)
    axis('off')
    imshow(im)
    
    # Sobel derivative filters
    
    imx = zeros(im.shape)
    filters.sobel(im, 1, imx)
    subplot(1, 4, 2)
    axis('off')
    imshow(imx)
    
    imy = zeros(im.shape)
    filters.sobel(im, 0, imy)
    subplot(1, 4, 3)
    axis('off')
    imshow(imy)
    
    #mag = numpy.sqrt(imx**2 + imy**2)
    mag = 255-sqrt(imx**2 + imy**2)
    subplot(1, 4, 4)
    axis('off')
    imshow(mag)
    
    show()
    

    得到图像:
    分别是原始灰度图像,x导数图像,y导数图像,梯度大小图像
    上面的方法有些缺陷,滤波器的尺寸需要随着图像分辨率的变化而变化。为了在任意尺度计算导数,我们可以使用高斯导数滤波器:

    def imx(im, sigma):
     imgx = zeros(im.shape)
     filters.gaussian_filter(im, sigma, (0, 1), imgx)
     return imgx
    
    def imy(im, sigma):
     imgy = zeros(im.shape)
     filters.gaussian_filter(im, sigma, (1, 0), imgy)
     return imgy
    
    def mag(im, sigma):
    
    # there's also gaussian_gradient_magnitude()
    
     #mag = numpy.sqrt(imgx**2 + imgy**2)
     imgmag = 255 - sqrt(imgx ** 2 + imgy ** 2)
     return imgmag
    
    im = array(Image.open('../data/empire.jpg').convert('L'))
    figure()
    gray()
    
    sigma = [2, 5, 10]
    
    for i in sigma:
     subplot(3, 4, 4*(sigma.index(i))+1)
     axis('off')
     imshow(im)
     imgx=imx(im, i)
     subplot(3, 4, 4*(sigma.index(i))+2)
     axis('off')
     imshow(imgx)
     imgy=imy(im, i)
     subplot(3, 4, 4*(sigma.index(i))+3)
     axis('off')
     imshow(imgy)
     imgmag=mag(im, i)
     subplot(3, 4, 4*(sigma.index(i))+4)
     axis('off')
     imshow(imgmag)
    
    show()
    

    得到图像:

    上中下分别为x导数、y导数、梯度大小图像。第一个为原始图像,二三四分别为σ等于2,5,10的高斯导数滤波器处理过的图像。

  3. 形态学:对象计数

    形态学通常用于处理二值图像,也能用于灰度图像。

    二值图像指图像的每个像素只能取0或1.

    二值图像通常是在计算物体的数目或者度量其大小时,对一副图像进行阈值化后的结果。

    figure()
    gray()
    im = array(Image.open('../data/houses.png').convert('L'))
    subplot(221)
    imshow(im)
    axis('off')
    
    #对图像进行阈值化操作,处理为二进制图像
    im = 1*(im < 128)
    
    #调用label方法统计个数
    labels, nbr_objects = measurements.label(im)
    print "Number of objects:", nbr_objects
    subplot(222)
    imshow(labels)
    axis('off')
    
    #由于原图有些对象间有微小连接,使用开操作去除。
    im_open = morphology.binary_opening(im, ones((9, 5)), iterations=2)
    subplot(223)
    imshow(im_open)
    axis('off')
    
    labels_open, nbr_objects_open = measurements.label(im_open)
    print "Number of objects:", nbr_objects_open
    subplot(224)
    imshow(labels_open)
    axis('off')
    
    show()
    

    使用开操作前后的统计个数分别为:

    Number of objects: 45
    Number of objects: 48
    

    图像:
    上面两图为开操作之前的,右边为label图。下图把那些微小连接去掉了。

  4. 用mat或图片保存数组。

    #打开mat文件
    data = scipy.io.loadmat('test.mat')
    ​``````language-text
    #创建字典
    data = {}
    #将变量x保存在字典中
    data['x'] = x
    scipy.io.savemat('test.mat',data)
    
    #将数组im保存到图片
    from scipy.misc import imsave
    imsave('test.jpg',im)
    
  5. 图像降噪。使用ROF模型。

    #原图
    im = array(Image.open('../data/empire.jpg').convert('L'))
    #使用ROF模型去噪
    U,T = rof.denoise(im,im)
    #使用高斯模糊
    G = filters.gaussian_filter(im,10)
    
    # save the result
    
    #imsave('synth_original.pdf',im)
    #imsave('synth_rof.pdf',U)
    #imsave('synth_gaussian.pdf',G)
    
    # plot
    
    figure()
    gray()
    
    subplot(1,3,1)
    imshow(im)
    #axis('equal')
    axis('off')
    
    subplot(1,3,2)
    imshow(G)
    #axis('equal')
    axis('off')
    
    subplot(1,3,3)
    imshow(U)
    #axis('equal')
    axis('off')
    
    show()
    

    得到图片:

    ROF去噪使得图像更平滑但是又结构清晰。可以对比高斯模糊去噪的效果。

局部图像描述子

介绍两种局部描述子算法。Harris角点检测算法和SIFT特征.

Harris角点检测器

如果像素周围显示存在多于一个方向的边,就认为该点为兴趣点。该点称为角点。

from pylab import *
from PIL import Image
from PCV.localdescriptors import harris

"""
Example of detecting Harris corner points (Figure 2-1 in the book).
"""

# 读入图像
im = array(Image.open('../data/empire.jpg').convert('L'))

# 检测harris角点
harrisim = harris.compute_harris_response(im)

# Harris响应函数
harrisim1 = 255 - harrisim

figure()
gray()

#画出Harris响应图
subplot(141)
imshow(harrisim1)
print harrisim1.shape
axis('off')
axis('equal')

threshold = [0.01, 0.05, 0.1]
for i, thres in enumerate(threshold):
    filtered_coords = harris.get_harris_points(harrisim, 6, thres)
    subplot(1, 4, i+2)
    imshow(im)
    print im.shape
    plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
    axis('off')

#原书采用的PCV中PCV harris模块
#harris.plot_harris_points(im, filtered_coords)

# plot only 200 strongest
# harris.plot_harris_points(im, filtered_coords[:200])

show()

如上,首先调用compute_harris_response检测Harris角点,它内部调用了高斯滤波器来抑制噪声强度。它返回像素值为Harris响应函数值的一幅图像。
然后我们从这幅图中挑选所需信息。在harris.get_harris_points中,我们选取像素值高于阈值的所有图像点,再加上额外限制,即角点之间的距离必须大于设定的最小距离。
然后将返回的这些点标注在原图上。
运行结果:
第一张图为Harris响应函数,后面的分别为阈值0.01,0.05,0.1检测出的角点。

在图像间寻找对应点

Harris角点检测器仅仅能检测图像中的兴趣点,但是没有给出通过比较图像间的兴趣点来寻找匹配角点的方法。我们需要在每个点加上描述子信息。
兴趣点描述子是分配给兴趣点的一个向量,描述该点附近的图像的表观信息。描述子越好,寻找到的对应点越好。我们用对应点或者点的对应来描述相同物体和场景点在不同图像上形成的像素点。
Harris角点的描述子通常是由周围图像像素块的灰度值,以及用于比较的归一化互相关矩阵构成的。图像的像素块由该像素点为中心的周围矩形部分图像构成。
该章代码在Harris中,首先我们调用get_descriptors方法来获取图像中的每个角点描述子。该函数将图像块像素压平为一个向量添加到描述子列表中。
然后调用match函数,来对应第一幅图中的每个角点描述子,使用归一化互相关,选取他们在第二幅图像中的匹配角点。
这些匹配可以通过在两边分别绘制图像,使用线段连接匹配的像素点来直观可视化。调用函数appendimages和plot_matches来完成。

# Figure 2-2下面的图

im1 = array(Image.open("../data/sf_view1.jpg").convert("L"))
im2 = array(Image.open("../data/sf_view2.jpg").convert("L"))

# resize加快匹配速度

im1 = imresize(im1, (im1.shape[1]/2, im1.shape[0]/2))
im2 = imresize(im2, (im2.shape[1]/2, im2.shape[0]/2))

wid = 5
harrisim = harris.compute_harris_response(im1, 5)
filtered_coords1 = harris.get_harris_points(harrisim, wid+1)
d1 = harris.get_descriptors(im1, filtered_coords1, wid)

harrisim = harris.compute_harris_response(im2, 5)
filtered_coords2 = harris.get_harris_points(harrisim, wid+1)
d2 = harris.get_descriptors(im2, filtered_coords2, wid)

print 'starting matching'
matches = harris.match_twosided(d1, d2)

figure()
gray()
harris.plot_matches(im1, im2, filtered_coords1, filtered_coords2, matches)
show()

得到图像:

从上图可以看到,存在一些不正确的匹配。而且不具有尺度不变性和旋转不变性,算法中像素块的大小也会影响对应匹配的结果。

SIFT 尺度不变特征变换

SIFT是近来最成功的图像局部描述子之一。SIFT描述子具有非常强的稳健性,它对于尺度、旋转和亮度都具有不变性,因此它可以用于三维视角和噪声的可靠匹配。

SIFT特征使用高斯差分函数来定位兴趣点。
为了实现旋转不变性,基于每个点周围图像梯度的方向和大小,SIFT描述子又引入了参考方向。SIFT描述子使用主方向描述参考方向。主方向使用方向直方图来度量。
为了对图像亮度具有稳健性,使用了图像梯度。在每个像素点附近选取子区域网格,在每个子区域内计算图像梯度方向直方图。每个子区域的直方图凭借起来组成描述子向量。
我们使用VLFeat提供的二进制文件来计算图像的SIFT特征。
因为SIFT需要图像格式为灰度.pgm来转换为sift文件,所以我们首先使用process_image函数来把图片转为灰度.pgm文件,调用sift命令转为目标sift文件。
sift文件文本类似如下:

274.094 177.346 0.973769 3.81914 18 2 0 35 138 10 6 35 76 4 0 0 1 0 20 138 114 0 0 0 0 0 8 112 11 0 0 0 0 0 1 7 14 11 5 77 138 25 5 6 83 16 2 9 5 14 10 66 138 9 0 0 0 1 2 119 41 2 0 0 0 0 0 17 14 1 4 79 50 40 22 82 104 26 1 17 41 48 6 7 138 83 0 0 0 1 0 9 43 11 0 0 0 0 0 4 20 1 0 0 5 48 51 44 138 48 0 0 5 37 8 8 138 49 0 0 0 0 0 14 24 1 0 0 0 0 0 5
290.517 177.408 0.923653 5.62849 8 10 8 3 11 10 5 13 34 15 1 2 9 10 3 22 159 111 0 0 0 0 0 72 125 33 0 0 0 2 1 33 30 3 5 8 11 7 11 38 115 20 15 35 5 4 0 12 159 136 0 1 0 0 0 32 89 50 0 0 0 0 1 8 57 7 13 12 11 1 1 22 116 15 6 26 7 7 5 24 159 59 0 2 1 1 0 16 80 19 0 0 0 0 1 5 13 2 6 4 2 11 6 13 117 64 0 0 0 6 5 26 159 150 0 0 0 0 0 0 28 15 0 0 0 0 0 0
......

每行前四个数据一次为兴趣点坐标、尺度和方向角度,后面的是对应描述符的128维向量。
然后read_features_from_file从SIFT文件读取特征属性值,将位置特征,描述子以矩阵形式返回。
最后使用plot_features函数来将特征点描绘到图片上。

imname = '../data/empire.jpg'
im = array(Image.open(imname).convert('L'))
sift.process_image(imname, 'empire.sift')
l1, d1 = sift.read_features_from_file('empire.sift')

figure()
gray()
subplot(131)
sift.plot_features(im, l1, circle=False)
subplot(132)
sift.plot_features(im, l1, circle=True)

# 检测harris角点

harrisim = harris.compute_harris_response(im)

subplot(133)
filtered_coords = harris.get_harris_points(harrisim, 6, 0.1)
imshow(im)
plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*')
axis('off')
show()

得到图像:

前面两个图是SIFT兴趣点,后面是Harris兴趣点。可以看到两个算法的点不同。

匹配描述子

将一幅图像中的特征匹配到另一幅图像的特征,一种稳健的准则是使用这两个特征距离和两个最匹配特征距离的比率。通过match函数,对第一幅图像中的每个描述子选取其在第二幅图像中的匹配。
为了增加稳健,我们会反过来再执行一次,最后我们仅保留满足这两种匹配准则的对应,使用match_twosided函数来实现。
然后和上面的Harris一样,使用appendimages和plot_matches函数来描绘出点。

if len(sys.argv) >= 3:
 im1f, im2f = sys.argv[1], sys.argv[2]
else:

# im1f = '../data/sf_view1.jpg'

# im2f = '../data/sf_view2.jpg'

 im1f = '../data/crans_1_small.jpg'
 im2f = '../data/crans_2_small.jpg'

# im1f = '../data/climbing_1_small.jpg'

# im2f = '../data/climbing_2_small.jpg'

im1 = array(Image.open(im1f))
im2 = array(Image.open(im2f))

sift.process_image(im1f, 'out_sift_1.txt')
l1, d1 = sift.read_features_from_file('out_sift_1.txt')
figure()
gray()
subplot(121)
sift.plot_features(im1, l1, circle=False)

sift.process_image(im2f, 'out_sift_2.txt')
l2, d2 = sift.read_features_from_file('out_sift_2.txt')
subplot(122)
sift.plot_features(im2, l2, circle=False)

#matches = sift.match(d1, d2)
matches = sift.match_twosided(d1, d2)
print '{} matches'.format(len(matches.nonzero()[0]))

figure()
gray()
sift.plot_matches(im1, im2, l1, l2, matches, show_below=True)
show()

得到图像:

匹配地理标记图片

我从google下载了一批白宫图片,放到一个目录中。然后使用局部描述子匹配,就是依次嵌套循环两两比对,计算出两图的匹配数量。

会生成矩阵数组,每对图像间的匹配特征数保存在matchscores数组中,其实是个对角线对称的矩阵,这样也只是为了美观。代码如下:

download_path = "/Users/shenzhongjie/Git/PCV/test/siftimages/" # set this to the path where you downloaded the panoramio images
path = "/Users/shenzhongjie/Git/PCV/test/siftimages/" # path to save thumbnails (pydot needs the full system path)

# list of downloaded filenames

imlist = imtools.get_imlist(download_path)
nbr_images = len(imlist)

# 将所有图像转为SIFT文件

featlist = [imname[:-3] + 'sift' for imname in imlist]
for i, imname in enumerate(imlist):
 sift.process_image(imname, featlist[i])

matchscores = zeros((nbr_images, nbr_images))

for i in range(nbr_images):
 for j in range(i, nbr_images): 

# 计算两两文件寻找匹配描述子,计算匹配数量生成矩阵

 print 'comparing ', imlist[i], imlist[j]
 l1, d1 = sift.read_features_from_file(featlist[i])
 l2, d2 = sift.read_features_from_file(featlist[j])
 matches = sift.match_twosided(d1, d2)
 nbr_matches = sum(matches > 0)
 print 'number of matches = ', nbr_matches
 matchscores[i, j] = nbr_matches
print "The match scores is: \n", matchscores

# copy values

for i in range(nbr_images):
 for j in range(i + 1, nbr_images): # no need to copy diagonal
 matchscores[j, i] = matchscores[i, j]

#可视化

threshold = 2 # min number of matches needed to create link

g = pydot.Dot(graph_type='graph') # don't want the default directed graph

for i in range(nbr_images):
 for j in range(i + 1, nbr_images):
 if matchscores[i, j] > threshold:

# first image in pair

 im = Image.open(imlist[i])
 im.thumbnail((100, 100))
 filename = path + str(i) + '.png'
 im.save(filename) # need temporary files of the right size
 g.add_node(pydot.Node(str(i), fontcolor='transparent', shape='rectangle', image=filename))

# second image in pair

 im = Image.open(imlist[j])
 im.thumbnail((100, 100))
 filename = path + str(j) + '.png'
 im.save(filename) # need temporary files of the right size
 g.add_node(pydot.Node(str(j), fontcolor='transparent', shape='rectangle', image=filename))

 g.add_edge(pydot.Edge(str(i), str(j)))
g.write_png('whitehouse.png')

从上面可以看到,我们得到矩阵后还要借助GraphViz图形库来绘制更直观的图,pydot是其的Python接口。还得借助pyparsing。分别安装,brew install graphviz,pip install pyparsing,pip install pydot.
最后得到的图如下:

我图片没选好,只下载了12张图片,得到这个图片。

图像聚类

K-means聚类

K-means是一种将输入数据划分为k个簇的简单的聚类算法。K-means反复提炼初始评估的类中心,步骤如下:

  1. 以随机或猜测的方式初始化类中心Ui,i=1…k;
  2. 将每个数据点归并到离它最近的类中心所属的类Ci;
  3. 对所有属于该类的数据点求平均,将平均值作为新的类中心;
  4. 重复步骤2和3直到收敛;

该算法是启发式提炼算法,普遍适用,但是不能更保证得到最优的结果。最大的缺陷是必须预先设定聚类数k,如果选择不恰当则会导致聚类出来的结果很差。

SciPy聚类包

使用scipy.cluster.vq中的K-means的实现。测试代码:

# 随机选取正态分布中值生成二维数组,[[-0.72053462 -1.98166101][-3.63685111 -0.26316214]..]

class1 = 1.5 * randn(100, 2)
class2 = randn(100, 2) + array([5, 5])

# 合并

features = vstack((class1, class2))

# 调用K-means计算默认20次,选择方差最小的结果

centroids, variance = kmeans(features, 2)

#对每个数据点进行归类,得到code进行归类
code, distance = vq(features, centroids)
figure()
ndx = where(code == 0)[0]
plot(features[ndx, 0], features[ndx, 1], '*')
ndx = where(code == 1)[0]
plot(features[ndx, 0], features[ndx, 1], 'r.')

#绘制预测出来的类中心
plot(centroids[:, 0], centroids[:, 1], 'go')

#axis('off')
show()

得到图像:

基于PCA的图像聚类

利用PCA算主成分,利用前40个主成分进行投影,用投影系数作为每幅图的向量描述符。载入之前PCA生成的pkl文件,在主成分上对图像进行投影,然后用K-means进行聚类。

为了方便直接用,下面代码里把之前PCA的代码也加上了。

# Uses sparse pca codepath.

imlist = imtools.get_imlist('../data/a_thumbs/')

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

im = array(Image.open(imlist[0])) # open one image to get the size
m, n = im.shape[:2] # get the size of the images
imnbr = len(imlist) # get the number of images
print "The number of images is %d" % imnbr

# Create matrix to store all flattened images

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

# PCA降维

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

# 保存均值和主成分

#f = open('./a_pca_modes.pkl', 'wb')
f = open('./a_pca_modes.pkl', 'wb')
pickle.dump(immean,f)
pickle.dump(V,f)
f.close()

# get list of images

imlist = imtools.get_imlist('../data/a_thumbs/')
imnbr = len(imlist)

# load model file

with open('./a_pca_modes.pkl','rb') as f:
 immean = pickle.load(f)
 V = pickle.load(f)

# create matrix to store all flattened images

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

# 投影到前40个主成分上

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

# k-means,k=4,四个类

projected = whiten(projected)
centroids,distortion = kmeans(projected,4)
code,distance = vq(projected,centroids)

# plot clusters

for k in range(4):
 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()

得到的结果如下四个图片:

在主成分上可视化图像

为了便于观察上面的利用PCA如何聚类的,我们可以在一对主成分方向的坐标上可视化这些图像。需要使用ImageDraw模块。代码如下:

imlist = imtools.get_imlist('../data/a_thumbs/')
imnbr = len(imlist)

# Load images, run PCA.

immatrix = array([array(Image.open(im)).flatten() for im in imlist], 'f')
V, S, immean = pca.pca(immatrix)

# Project on 2 PCs.

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

# height and width

h, w = 1200, 1200

# create a new image with a white background

img = Image.new('RGB', (w, h), (255, 255, 255))
draw = ImageDraw.Draw(img)

# draw axis

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

# scale coordinates to fit

scale = abs(projected).max(0)
scaled = floor(array([(p/scale) * (w/2 - 20, h/2 - 20) + (w/2, h/2)
 for p in projected])).astype(int)

# paste thumbnail of each image

for i in range(imnbr):
 nodeim = Image.open(imlist[i])
 nodeim.thumbnail((25, 25))
 ns = nodeim.size
 box = (scaled[i][0] - ns[0] // 2, scaled[i][1] - ns[1] // 2,
 scaled[i][0] + ns[0] // 2 + 1, scaled[i][1] + ns[1] // 2 + 1)
 img.paste(nodeim, box)

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

figure()
imshow(img)
axis('off')
img.save('pca_font.png')
show()

我得到的图片如下:
坐标的左上角为(0,0),两个图距离越近相似度越高。

像素聚类

对单幅图像中的像素而非全部图像进行聚类的例子,将图像区域或像素合并成有意义的部分称为_图像分割_。现在我们只谈论在RGB三通道的像素值上运用K-means进行聚类,分割暂不说。
代码用一个步长为steps的方形网格在图像中滑动,每滑一次对网格中图像区域像素求平均值,将其作为新生成的低分辨率图像对应位置处的像素值,并用K-means进行聚类。
K-means的输入时一个有steps*steps行的数组,数组的每一行有三列,每列分别为区域块RGB三个通道的像素平均值。

def clusterpixels(infile, k, steps):
 im = array(Image.open(infile))
 dx = im.shape[0] / steps
 dy = im.shape[1] / steps

# compute color features for each region

 features = []
 for x in range(steps):
 for y in range(steps):
 R = mean(im[x * dx:(x + 1) * dx, y * dy:(y + 1) * dy, 0])
 G = mean(im[x * dx:(x + 1) * dx, y * dy:(y + 1) * dy, 1])
 B = mean(im[x * dx:(x + 1) * dx, y * dy:(y + 1) * dy, 2])
 features.append([R, G, B])
 features = array(features, 'f') # make into array

# 聚类, k是聚类数目

 centroids, variance = kmeans(features, k)
 code, distance = vq(features, centroids)

# create image with cluster labels

 codeim = code.reshape(steps, steps)
 codeim = imresize(codeim, im.shape[:2])
 return codeim

k=3
infile_empire = '../data/empire.jpg'
im_empire = array(Image.open(infile_empire))
infile_boy_on_hill = '../data/boy_on_hill.jpg'
im_boy_on_hill = array(Image.open(infile_boy_on_hill))
steps = (50, 100) # image is divided in steps*steps region
print steps[0], steps[-1]

#显示原图empire.jpg
figure()
subplot(231)
axis('off')
imshow(im_empire)

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

codeim= clusterpixels(infile_empire, k, steps[0])
subplot(232)
#ax1.set_title('Image')
axis('off')
imshow(codeim)

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

codeim= clusterpixels(infile_empire, k, steps[-1])
ax1 = subplot(233)
#ax1.set_title('Image')
axis('off')
imshow(codeim)

#显示原图empire.jpg
subplot(234)
axis('off')
imshow(im_boy_on_hill)

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

codeim= clusterpixels(infile_boy_on_hill, k, steps[0])
subplot(235)
#ax1.set_title('Image')
axis('off')
imshow(codeim)

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

codeim= clusterpixels(infile_boy_on_hill, k, steps[-1])
subplot(236)
axis('off')
imshow(codeim)

show()

得到图片:
这儿变形了,因为调用resize时报了一个参数个数错误,所以去了一个参数。

基于层次聚类的图像聚类

层次聚类(凝聚式聚类)是一种简单有效的聚类算法,其思想是基于样本间成对距离建立一个简相似性树。该算法首先将特征向量距离最近的两个样本归并为一组,并在树中创建一个“平均”节点,将这两个距离最近的样本作为该“平均”节点下的子节点;然后在剩下的包含任意平均节点的样本中寻找下一个最近的树,重复进行前面的操作,在每一个节点处保存了两个子节点之间的距离。遍历整个树,通过设定的阈值,遍历过程可以在比阈值大的节点位置终止,从而提取出聚类簇。
层次聚类有若干优点,如利用树结构可以可可视化数据间的关系,并显示这些簇是如何关联的。在树中,一个好的特征向量可以给出一个很好的分离结果。另外一个优点是,对于给定的不同的阈值,可以直接利用原来的树。不足时我们需要选择一个合适的阈值。

类似上面的K-means,我们先创建一些二位数据点,如下进行聚类:

class1 = 1.5 * randn(100,2)
class2 = randn(100,2) + array([5,5])
features = vstack((class1,class2))

tree = hcluster.hcluster(features)
clusters = tree.extract_clusters(5)
print 'number of clusters', len(clusters)
for c in clusters:
 print c.get_cluster_elements()

得到两个聚类簇,按理来说一个大于100,一个小于100。不过我得到的有点误差:

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

现在来试试引用于图像聚类。这儿我们用颜色直方图作为每幅图像的特征向量。
我们将RGB三个颜色通道作为特征向量,将其传递到NumPy的histogramdd()中,该函数能够计算多维直方图。
我们在每个颜色通道中使用8个小区间进行量化,将三个通道量化后的小区间拉成一行后便可用512(888)维的特征向量描述每幅图像。为了避免图像尺寸不一致我们用normed=True归一化直方图,并将每个颜色通道范围设置为0…255.将reshape()第一个参数设置为-1会自动确定正确地尺寸,故可以创建一个输入数组来计算以RGB颜色值为行向量的直方图。
代码如下:

# create a list of images

path = '../data/flickr-sunsets-small/'
imlist = [os.path.join(path, f) for f in os.listdir(path) if f.endswith('.jpg')]

# extract feature vector (8 bins per color channel)

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

# multi-dimensional histogram

 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)

# visualize clusters with some (arbitrary) threshold

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

# plot images for clusters with more than 3 elements

for c in clusters:
 elements = c.get_cluster_elements()
 nbr_elements = len(elements)
 if nbr_elements > 3:
 figure()
 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='sunset.pdf')

得到结果:

谱聚类

谱聚类和K-means和层次聚类完全不同。
对于n个元素(如n幅图像),相似矩阵是一个n*n的矩阵,矩阵每个元素表示两两之间的相似性分数。谱聚类是有相似性矩阵构建谱矩阵而得名的。对该谱矩阵进行特征向量可以用于降维,然后聚类。
其优点之一是仅需输入相似性矩阵,并且可以采用你所想到的任意度量方式构建该相似性矩阵。像K-means和层次聚类需要计算那些特征向量求平均;为了计算平均值,会将特征或描述子限制为向量,而对于谱方法,特征向量就没有类别限制,只要有一个距离或相似性的概念就可以。

imlist = imtools.get_imlist('../data/a_thumbs')
imnbr = len(imlist)

# Load images, run PCA.

immatrix = array([array(Image.open(im)).flatten() for im in imlist], 'f')
V, S, immean = pca.pca(immatrix)

# Project on 2 PCs.

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

n = len(projected)

# 计算距离矩阵

S = array([[ sqrt(sum((projected[i]-projected[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))

# 计算矩阵L的特征向量

U,sigma,V = linalg.svd(L)
k = 5

# 从矩阵L的前k个特征向量中创建特征向量

# 叠加特征向量作为数组的列

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

# k-means

features = whiten(features)
centroids,distortion = kmeans(features,k)
code,distance = vq(features,centroids)

# plot clusters

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

得到如下聚类结果:

上面的结果是用拉普拉斯矩阵的特征向量对字体图像进行谱聚类。

我们也可以在没有任何特征向量或是没有严格定义相似性的例子中尝试该算法。比如之前的白宫例子中我们用多少匹配的描述子来计算的。前面我们在匹配描述子中生成过一个用分数表示的相似性矩阵,这儿可以类似转换一下进行谱聚类。

# list of downloaded filenames

imlist = imtools.get_imlist('./siftimages')
nbr_images = len(imlist)

# extract features

#featlist = [imname[:-3] + 'sift' for imname in imlist]
#for i, imname in enumerate(imlist):

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

featlist = glob.glob('./siftimages/*.sift')

matchscores = zeros((nbr_images, nbr_images))

for i in range(nbr_images):
 for j in range(i, nbr_images): # only compute upper triangle
 print 'comparing ', imlist[i], imlist[j]
 l1, d1 = sift.read_features_from_file(featlist[i])
 l2, d2 = sift.read_features_from_file(featlist[j])
 matches = sift.match_twosided(d1, d2)
 nbr_matches = sum(matches > 0)
 print 'number of matches = ', nbr_matches
 matchscores[i, j] = nbr_matches
print "The match scores is: \n", matchscores

# copy values

for i in range(nbr_images):
 for j in range(i + 1, nbr_images): # no need to copy diagonal
 matchscores[j, i] = matchscores[i, j]

n = len(imlist)

# load the similarity matrix and reformat

S = matchscores
S = 1 / (S + 1e-6)

# create Laplacian matrix

rowsum = sum(S,axis=0)
D = diag(1 / sqrt(rowsum))
I = identity(n)
L = I - dot(D,dot(S,D))

# compute eigenvectors of L

U,sigma,V = linalg.svd(L)
k = 10

# create feature vector from k first eigenvectors

# by stacking eigenvectors as columns

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

# k-means

features = whiten(features)
centroids,distortion = kmeans(features,k)
code,distance = vq(features,centroids)

# plot clusters

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

这是k=2的情况:

不是特别准确。看看k=10的情况:

还是比较糟糕,大多是单个的,去除后的也不太理想。

图像搜索

在大型图像数据库上,CBIR技术用于检索在视觉上具有相似性的图像。这样返回的图像可以是颜色相似、纹理相似、图像中的物体或场景相似。

从文本挖掘中获取灵感——矢量空间模型

矢量空间模型是一个用于表示和搜索文本文档的模型。名字来源于用矢量表示文本文档,这些矢量是由文本词频直方图构成的。矢量包含了每个单词出现的次数,而且在其他别的地方包含很多0元素。由于其忽略了单词出现的顺序和位置,该模型也被称为BOW来表示模型。

通过单词计数来构建文档直方图向量v,从而建立文档索引。
会忽略一些词,比如“这”,“是”,这些词被称为停用词。
由于每篇文档长度不同,故除以直方图总将向量归一化为单位长度。
对于直方图向量中的每个元素,一般根据每个单词的重要性来赋予权重。
一般数据集(语料库)中一个单词的重要性和他在文档中出现的次数成正比,与他在语料库中出现的次数成反比。

最常用的权重是tf-idf(词频-逆向文档词频)。

视觉单词

为了将文本挖掘技术应用到图像中,我们需要建立视觉等效单词;
使用SIFT局部描述子,它的思想是将描述子空间量化成一些典型实例,并将图像中的每个描述子指派到其他的某个实例中。
这些实例可以通过分析训练图像集确定,被称为视觉单词。
所有这些视觉单词构成的集合称为视觉词汇或视觉码本。

从一个很大的训练图像集提取特征描述子,利用一些聚类算法可以构建出视觉单词。通常用K-means。
视觉单词只是在给定特征描述子空间中的一组向量集。
用视觉单词直方图来表示图像,则该模型便成为BOW模型。