转自 https://my.oschina.net/ranvane/blog/352341

 

Image manipulation and processing using Numpy and Scipy

翻译自：http://scipy-lectures.github.com/advanced/image_processing/index.html

作者：Emmanuelle Gouillart, Gaël Varoquaux

1. 图像 = 2-D 数值数组

2.  

3. (或者 3-D: CT, MRI, 2D + 时间; 4-D, ...)

4.  

5. 这里 图像 == Numpy数组 np.array

这个教程中使用的工具：

• numpy：基本数组操作

• scipy：scipy.ndimage子模块致力于图像处理(n维图像)。参见http://docs.scipy.org/doc/scipy/reference/tutorial/ndimage.html

  from scipy import ndimage

• 一些例子用到了使用np.array的特殊的工具箱：

  • Scikit Image

  • scikit-learn

图像中的常见问题有：

• 输入/输出，呈现图像

• 基本操作：裁剪、翻转、旋转……

• 图像滤镜：消噪，锐化

• 图像分割：不同对应对象的像素标记

更有力和完整的模块：

• OpenCV (Python绑定)

• CellProfiler

• ITK,Python绑定

• 更多……

目录

• 打开和读写图像文件

• 呈现图像

• 基本操作

  • 统计信息

  • 几何转换

• 图像滤镜

  • 模糊/平滑

  • 锐化

  • 消噪

  • 数学形态学

• 特征提取

  • 边缘检测

  • 分割

• 测量对象属性：ndimage.measurements

• Footnotes

打开和读写图像文件

将一个数组写入文件：

1. In [1]: from scipy import misc

2.  

3. In [2]: l = misc.lena()

4.  

5. In [3]: misc.imsave('lena.png', l)  # uses the Image module (PIL)

6.  

7. In [4]: import pylab as pl

8.  

9. In [5]: pl.imshow(l)

10. Out[5]: <matplotlib.image.AxesImage at 0x4118110>

从一个图像文件创建数组：

1. In [7]: lena = misc.imread('lena.png')

2.  

3. In [8]: type(lena)

4. Out[8]: numpy.ndarray

5.  

6. In [9]: lena.shape, lena.dtype

7. Out[9]: ((512, 512), dtype('uint8'))

8位图像(0-255)的dtype是uint8

打开一个raw文件(相机， 3-D图像)

1. In [10]: l.tofile('lena.raw')  # 创建一个raw文件

2.  

3. In [14]: lena_from_raw = np.fromfile('lena.raw', dtype=np.int64)

4.  

5. In [15]: lena_from_raw.shape

6. Out[15]: (262144,)

7.  

8. In [16]: lena_from_raw.shape = (512, 512)

9.  

10. In [17]: import os

11.  

12. In [18]: os.remove('lena.raw')

需要知道图像的shape和dtype(如何区分隔数据字节)

对于大数据，使用np.memmap进行内存映射：

In [21]: lena_memmap = np.memmap('lena.raw', dtype=np.int64, shape=(512,512))

(数据从文件读取，而不是载入内存)

处理一个列表的图像文件：

1. In [22]: for i in range(10):

2.    ....:     im = np.random.random_integers(0, 255, 10000).reshape((100, 100))

3.    ....:     misc.imsave('random_%02d.png' % i, im)

4.    ....:     

5.  

6. In [23]: from glob import glob

7.  

8. In [24]: filelist = glob('random*.png')

9.  

10. In [25]: filelist.sort()

呈现图像

使用matplotlib和imshow将图像呈现在matplotlib图像(figure)中：

1. In [29]: l = misc.lena()

2.  

3. In [30]: import matplotlib.pyplot as plt

4.  

5. In [31]: plt.imshow(l, cmap=plt.cm.gray)

6. Out[31]: <matplotlib.image.AxesImage at 0x4964990>

通过设置最大最小之增加对比：

1. In [33]: plt.imshow(l, cmap=plt.cm.gray, vmin=30, vmax=200)

2. Out[33]: <matplotlib.image.AxesImage at 0x50cb790>

3.  

4. In [34]: plt.axis('off')  # 移除axes和ticks

5. Out[34]: (-0.5, 511.5, 511.5, -0.5)

绘制等高线：1

ln[7]: plt.contour(l, [60, 211])

更好地观察强度变化，使用interpolate=‘nearest’：

1. In [7]: plt.imshow(l[200:220, 200:220], cmap=plt.cm.gray)

2. Out[7]: <matplotlib.image.AxesImage at 0x3bbe610>

3.  

4. In [8]: plt.imshow(l[200:220, 200:220], cmap=plt.cm.gray, interpolation='nearest')

5. Out[8]: <matplotlib.image.AxesImage at 0x3ed3250>

其它包有时使用图形工具箱来可视化(GTK，Qt)：2

1. In [9]: import skimage.io as im_io

2.  

3. In [21]: im_io.use_plugin('gtk', 'imshow')

4.  

5. In [22]: im_io.imshow(l)

3-D可视化：Mayavi

参见可用Mayavi进行3-D绘图和体积数据

• 图形平面工具

• 等值面

• ……

基本操作

图像是数组：使用整个numpy机理。

1. >>> lena = misc.lena()

2. >>> lena[0, 40]

3. 166

4. >>> # Slicing

5. >>> lena[10:13, 20:23]

6. array([[158, 156, 157],

7. [157, 155, 155],

8. [157, 157, 158]])

9. >>> lena[100:120] = 255

10. >>>

11. >>> lx, ly = lena.shape

12. >>> X, Y = np.ogrid[0:lx, 0:ly]

13. >>> mask = (X - lx/2)**2 + (Y - ly/2)**2 > lx*ly/4

14. >>> # Masks

15. >>> lena[mask] = 0

16. >>> # Fancy indexing

17. >>> lena[range(400), range(400)] = 255

统计信息

1. >>> lena = scipy.lena()

2. >>> lena.mean()

3. 124.04678344726562

4. >>> lena.max(), lena.min()

5. (245, 25)

np.histogram

几何转换

1. >>> lena = scipy.lena()

2. >>> lx, ly = lena.shape

3. >>> # Cropping

4. >>> crop_lena = lena[lx/4:-lx/4, ly/4:-ly/4]

5. >>> # up <-> down flip

6. >>> flip_ud_lena = np.flipud(lena)

7. >>> # rotation

8. >>> rotate_lena = ndimage.rotate(lena, 45)

9. >>> rotate_lena_noreshape = ndimage.rotate(lena, 45, reshape=False)

示例源码

图像滤镜

局部滤镜：用相邻像素值的函数替代当前像素的值。

相邻：方形(指定大小)，圆形， 或者更多复杂的_结构元素_。

模糊/平滑

scipy.ndimage中的_高斯滤镜_：

1. >>> from scipy import misc

2. >>> from scipy import ndimage

3. >>> lena = misc.lena()

4. >>> blurred_lena = ndimage.gaussian_filter(lena, sigma=3)

5. >>> very_blurred = ndimage.gaussian_filter(lena, sigma=5)

均匀滤镜

>>> local_mean = ndimage.uniform_filter(lena, size=11)

示例源码

锐化

锐化模糊图像：

1. >>> from scipy import misc

2. >>> lena = misc.lena()

3. >>> blurred_l = ndimage.gaussian_filter(lena, 3)

通过增加拉普拉斯近似增加边缘权重：

1. >>> filter_blurred_l = ndimage.gaussian_filter(blurred_l, 1)

2. >>> alpha = 30

3. >>> sharpened = blurred_l + alpha * (blurred_l - filter_blurred_l)

示例源码

消噪

向lena增加噪声：

1. >>> from scipy import misc

2. >>> l = misc.lena()

3. >>> l = l[230:310, 210:350]

4. >>> noisy = l + 0.4*l.std()*np.random.random(l.shape)

_高斯滤镜_平滑掉噪声……还有边缘：

>>> gauss_denoised = ndimage.gaussian_filter(noisy, 2)

大多局部线性各向同性滤镜都模糊图像(ndimage.uniform_filter)

_中值滤镜_更好地保留边缘：

>>> med_denoised = ndimage.median_filter(noisy, 3)

示例源码

中值滤镜：对直边界效果更好(低曲率)：

1. >>> im = np.zeros((20, 20))

2. >>> im[5:-5, 5:-5] = 1

3. >>> im = ndimage.distance_transform_bf(im)

4. >>> im_noise = im + 0.2*np.random.randn(*im.shape)

5. >>> im_med = ndimage.median_filter(im_noise, 3)

示例源码

其它排序滤波器：ndimage.maximum_filter,ndimage.percentile_filter

其它局部非线性滤波器：维纳滤波器(scipy.signal.wiener)等

非局部滤波器

_总变差(TV)_消噪。找到新的图像让图像的总变差(正态L1梯度的积分)变得最小，当接近测量图像时：

1. >>> # from skimage.filter import tv_denoise

2. >>> from tv_denoise import tv_denoise

3. >>> tv_denoised = tv_denoise(noisy, weight=10)

4. >>> # More denoising (to the expense of fidelity to data)

5. >>> tv_denoised = tv_denoise(noisy, weight=50)

总变差滤镜tv_denoise可以从skimage中获得，(文档:http://scikit-image.org/docs/dev/api/skimage.filter.html#denoise-tv)，但是为了方便我们在这个教程中作为一个_单独模块_导入。

示例源码

数学形态学

参见：http://en.wikipedia.org/wiki/Mathematical_morphology

结构元素：

1. >>> el = ndimage.generate_binary_structure(2, 1)

2. >>> el

3. array([[False,  True, False],

4.        [ True,  True,  True],

5.        [False,  True, False]], dtype=bool)

6. >>> el.astype(np.int)

7. array([[0, 1, 0],

8.        [1, 1, 1],

9.        [0, 1, 0]])

腐蚀 = 最小化滤镜。用结构元素覆盖的像素的最小值替代一个像素值：

1. >>> a = np.zeros((7,7), dtype=np.int)

2. >>> a[1:6, 2:5] = 1

3. >>> a

4. array([[0, 0, 0, 0, 0, 0, 0],

5.        [0, 0, 1, 1, 1, 0, 0],

6.        [0, 0, 1, 1, 1, 0, 0],

7.        [0, 0, 1, 1, 1, 0, 0],

8.        [0, 0, 1, 1, 1, 0, 0],

9.        [0, 0, 1, 1, 1, 0, 0],

10.        [0, 0, 0, 0, 0, 0, 0]])

11. >>> ndimage.binary_erosion(a).astype(a.dtype)

12. array([[0, 0, 0, 0, 0, 0, 0],

13.        [0, 0, 0, 0, 0, 0, 0],

14.        [0, 0, 0, 1, 0, 0, 0],

15.        [0, 0, 0, 1, 0, 0, 0],

16.        [0, 0, 0, 1, 0, 0, 0],

17.        [0, 0, 0, 0, 0, 0, 0],

18.        [0, 0, 0, 0, 0, 0, 0]])

19. >>> #Erosion removes objects smaller than the structure

20. >>> ndimage.binary_erosion(a, structure=np.ones((5,5))).astype(a.dtype)

21. array([[0, 0, 0, 0, 0, 0, 0],

22.        [0, 0, 0, 0, 0, 0, 0],

23.        [0, 0, 0, 0, 0, 0, 0],

24.        [0, 0, 0, 0, 0, 0, 0],

25.        [0, 0, 0, 0, 0, 0, 0],

26.        [0, 0, 0, 0, 0, 0, 0],

27.        [0, 0, 0, 0, 0, 0, 0]])

膨胀：最大化滤镜:

1. >>> a = np.zeros((5, 5))

2. >>> a[2, 2] = 1

3. >>> a

4. array([[ 0.,  0.,  0.,  0.,  0.],

5.        [ 0.,  0.,  0.,  0.,  0.],

6.        [ 0.,  0.,  1.,  0.,  0.],

7.        [ 0.,  0.,  0.,  0.,  0.],

8.        [ 0.,  0.,  0.,  0.,  0.]])

9. >>> ndimage.binary_dilation(a).astype(a.dtype)

10. array([[ 0.,  0.,  0.,  0.,  0.],

11.        [ 0.,  0.,  1.,  0.,  0.],

12.        [ 0.,  1.,  1.,  1.,  0.],

13.        [ 0.,  0.,  1.,  0.,  0.],

14.        [ 0.,  0.,  0.,  0.,  0.]])

对灰度值图像也有效：

1. >>> np.random.seed(2)

2. >>> x, y = (63*np.random.random((2, 8))).astype(np.int)

3. >>> im[x, y] = np.arange(8)

4.  

5. >>> bigger_points = ndimage.grey_dilation(im, size=(5, 5), structure=np.ones((5, 5)))

6.  

7. >>> square = np.zeros((16, 16))

8. >>> square[4:-4, 4:-4] = 1

9. >>> dist = ndimage.distance_transform_bf(square)

10. >>> dilate_dist = ndimage.grey_dilation(dist, size=(3, 3), \

11. ...         structure=np.ones((3, 3)))

示例源码

开操作：腐蚀+膨胀：

应用：移除噪声

1. >>> square = np.zeros((32, 32))

2. >>> square[10:-10, 10:-10] = 1

3. >>> np.random.seed(2)

4. >>> x, y = (32*np.random.random((2, 20))).astype(np.int)

5. >>> square[x, y] = 1

6.  

7. >>> open_square = ndimage.binary_opening(square)

8.  

9. >>> eroded_square = ndimage.binary_erosion(square)

10. >>> reconstruction = ndimage.binary_propagation(eroded_square, mask=square)

示例源码

闭操作：膨胀+腐蚀

许多其它数学分形：击中(hit)和击不中(miss)变换，tophat等等。

特征提取

边缘检测

合成数据：

1. >>> im = np.zeros((256, 256))

2. >>> im[64:-64, 64:-64] = 1

3. >>>

4. >>> im = ndimage.rotate(im, 15, mode='constant')

5. >>> im = ndimage.gaussian_filter(im, 8)

使用_梯度操作(Sobel)_来找到搞强度的变化：

1. >>> sx = ndimage.sobel(im, axis=0, mode='constant')

2. >>> sy = ndimage.sobel(im, axis=1, mode='constant')

3. >>> sob = np.hypot(sx, sy)

示例源码

canny滤镜

Canny滤镜可以从skimage中获取(文档)，但是为了方便我们在这个教程中作为一个_单独模块_导入：

1. >>> #from skimage.filter import canny

2. >>> #or use module shipped with tutorial

3. >>> im += 0.1*np.random.random(im.shape)

4. >>> edges = canny(im, 1, 0.4, 0.2) # not enough smoothing

5. >>> edges = canny(im, 3, 0.3, 0.2) # better parameters

示例源码

需要调整几个参数……过度拟合的风险

分割

• 基于_直方图_的分割(没有空间信息)

  1.   >>> n = 10

  2.   >>> l = 256

  3.   >>> im = np.zeros((l, l))

  4.   >>> np.random.seed(1)

  5.   >>> points = l*np.random.random((2, n**2))

  6.   >>> im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1

  7.   >>> im = ndimage.gaussian_filter(im, sigma=l/(4.*n))

  8.     

  9.   >>> mask = (im > im.mean()).astype(np.float)

  10.   >>> mask += 0.1 * im

  11.   >>> img = mask + 0.2*np.random.randn(*mask.shape)

  12.     

  13.   >>> hist, bin_edges = np.histogram(img, bins=60)

  14.   >>> bin_centers = 0.5*(bin_edges[:-1] + bin_edges[1:])

  15.     

  16.   >>> binary_img = img > 0.5

示例源码

自动阈值：使用高斯混合模型：

1. >>> mask = (im > im.mean()).astype(np.float)

2. >>> mask += 0.1 * im

3. >>> img = mask + 0.3*np.random.randn(*mask.shape)

4.  

5. >>> from sklearn.mixture import GMM

6. >>> classif = GMM(n_components=2)

7. >>> classif.fit(img.reshape((img.size, 1))) 

8. GMM(...)

9.  

10. >>> classif.means_

11. array([[ 0.9353155 ],

12.        [-0.02966039]])

13. >>> np.sqrt(classif.covars_).ravel()

14. array([ 0.35074631,  0.28225327])

15. >>> classif.weights_

16. array([ 0.40989799,  0.59010201])

17. >>> threshold = np.mean(classif.means_)

18. >>> binary_img = img > threshold

使用数学形态学来清理结果：

1. >>> # Remove small white regions

2. >>> open_img = ndimage.binary_opening(binary_img)

3. >>> # Remove small black hole

4. >>> close_img = ndimage.binary_closing(open_img)

示例源码

练习

参看重建(reconstruction)操作(腐蚀+传播(propagation))产生比开/闭操作更好的结果：

1. >>> eroded_img = ndimage.binary_erosion(binary_img)

2. >>> reconstruct_img = ndimage.binary_propagation(eroded_img, mask=binary_img)

3. >>> tmp = np.logical_not(reconstruct_img)

4. >>> eroded_tmp = ndimage.binary_erosion(tmp)

5. >>> reconstruct_final = np.logical_not(ndimage.binary_propagation(eroded_tmp, mask=tmp))

6. >>> np.abs(mask - close_img).mean()

7. 0.014678955078125

8. >>> np.abs(mask - reconstruct_final).mean()

9. 0.0042572021484375

练习

检查首次消噪步骤(中值滤波，总变差)如何更改直方图，并且查看是否基于直方图的分割更加精准了。

• _基于图像_的分割：使用空间信息

  1.   >>> from sklearn.feature_extraction import image

  2.   >>> from sklearn.cluster import spectral_clustering

  3.     

  4.   >>> l = 100

  5.   >>> x, y = np.indices((l, l))

  6.     

  7.   >>> center1 = (28, 24)

  8.   >>> center2 = (40, 50)

  9.   >>> center3 = (67, 58)

  10.   >>> center4 = (24, 70)

  11.   >>> radius1, radius2, radius3, radius4 = 16, 14, 15, 14

  12.     

  13.   >>> circle1 = (x - center1[0])**2 + (y - center1[1])**2 < radius1**2

  14.   >>> circle2 = (x - center2[0])**2 + (y - center2[1])**2 < radius2**2

  15.   >>> circle3 = (x - center3[0])**2 + (y - center3[1])**2 < radius3**2

  16.   >>> circle4 = (x - center4[0])**2 + (y - center4[1])**2 < radius4**2

  17.     

  18.   >>> # 4 circles

  19.   >>> img = circle1 + circle2 + circle3 + circle4

  20.   >>> mask = img.astype(bool)

  21.   >>> img = img.astype(float)

  22.     

  23.   >>> img += 1 + 0.2*np.random.randn(*img.shape)

  24.   >>> # Convert the image into a graph with the value of the gradient on

  25.   >>> # the edges.

  26.   >>> graph = image.img_to_graph(img, mask=mask)

  27.     

  28.   >>> # Take a decreasing function of the gradient: we take it weakly

  29.   >>> # dependant from the gradient the segmentation is close to a voronoi

  30.   >>> graph.data = np.exp(-graph.data/graph.data.std())

  31.     

  32.   >>> labels = spectral_clustering(graph, k=4, mode='arpack')

  33.   >>> label_im = -np.ones(mask.shape)

  34.   >>> label_im[mask] = labels

---

测量对象属性：ndimage.measurements

合成数据：

1. >>> n = 10

2. >>> l = 256

3. >>> im = np.zeros((l, l))

4. >>> points = l*np.random.random((2, n**2))

5. >>> im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1

6. >>> im = ndimage.gaussian_filter(im, sigma=l/(4.*n))

7. >>> mask = im > im.mean()

• 连接成分分析

  标记连接成分：ndimage.label

  1.   >>> label_im, nb_labels = ndimage.label(mask)

  2.   >>> nb_labels # how many regions?

  3.   23

  4.   >>> plt.imshow(label_im)        

  5.   <matplotlib.image.AxesImage object at ...>

示例源码

计算每个区域的尺寸，均值等等：

1. >>> sizes = ndimage.sum(mask, label_im, range(nb_labels + 1))

2. >>> mean_vals = ndimage.sum(im, label_im, range(1, nb_labels + 1))

计算小的连接成分：

1. >>> mask_size = sizes < 1000

2. >>> remove_pixel = mask_size[label_im]

3. >>> remove_pixel.shape

4. (256, 256)

5. >>> label_im[remove_pixel] = 0

6. >>> plt.imshow(label_im)        

7. <matplotlib.image.AxesImage object at ...>

现在使用np.searchsorted重新分配标签：

1. >>> labels = np.unique(label_im)

2. >>> label_im = np.searchsorted(labels, label_im)

示例源码

找到关注的封闭对象区域：3

1. >>> slice_x, slice_y = ndimage.find_objects(label_im==4)[0]

2. >>> roi = im[slice_x, slice_y]

3. >>> plt.imshow(roi)     

4. <matplotlib.image.AxesImage object at ...>

示例源码

其它空间测量：ndiamge.center_of_mass,ndimage.maximum_position等等。

可以在分割应用限制范围之外使用。

示例：块平均(block mean)：

1. m scipy import misc

2. >>> l = misc.lena()

3. >>> sx, sy = l.shape

4. >>> X, Y = np.ogrid[0:sx, 0:sy]

5. >>> regions = sy/6 * (X/4) + Y/6  # note that we use broadcasting

6. >>> block_mean = ndimage.mean(l, labels=regions, index=np.arange(1,

7. ...     regions.max() +1))

8. >>> block_mean.shape = (sx/4, sy/6)

示例源码

当区域不是正则的4块状时，使用stride技巧更有效(示例：fake dimensions with strides)

非正则空间(Non-regular-spaced)区块：径向平均：

1. >>> sx, sy = l.shape

2. >>> X, Y = np.ogrid[0:sx, 0:sy]

3. >>> r = np.hypot(X - sx/2, Y - sy/2)

4. >>> rbin = (20* r/r.max()).astype(np.int)

5. >>> radial_mean = ndimage.mean(l, labels=rbin, index=np.arange(1, rbin.max() +1))

示例源码

• 其它测量

相关函数，傅里叶/小波谱等。

一个使用数学形态学的例子：粒度(http://en.wikipedia.org/wiki/Granulometry_%28morphology%29)

1. >>> def disk_structure(n):

2. ...     struct = np.zeros((2 * n + 1, 2 * n + 1))

3. ...     x, y = np.indices((2 * n + 1, 2 * n + 1))

4. ...     mask = (x - n)**2 + (y - n)**2 <= n**2

5. ...     struct[mask] = 1

6. ...     return struct.astype(np.bool)

7. ...

8. >>>

9. >>> def granulometry(data, sizes=None):

10. ...     s = max(data.shape)

11. ...     if sizes == None:

12. ...         sizes = range(1, s/2, 2)

13. ...     granulo = [ndimage.binary_opening(data, \

14. ...         structure=disk_structure(n)).sum() for n in sizes]

15. ...     return granulo

16. ...

17. >>>

18. >>> np.random.seed(1)

19. >>> n = 10

20. >>> l = 256

21. >>> im = np.zeros((l, l))

22. >>> points = l*np.random.random((2, n**2))

23. >>> im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1

24. >>> im = ndimage.gaussian_filter(im, sigma=l/(4.*n))

25. >>>

26. >>> mask = im > im.mean()

27. >>>

28. >>> granulo = granulometry(mask, sizes=np.arange(2, 19, 4))

示例源码

Footnotes

1. 占位 ↩

2. ValueError: can not convert int64 to uint8. ↩

3. 根据以上操作剩下的区域选择区域，因为是随机生成可能结果不通，label_im==4未必留下来了。 ↩

4. 正则空间 ↩

原文地址：http://reverland.org/python/2012/08/24/scipy/