2.6. Image manipulation and processing using Numpy and Scipy

作者Emmanuelle Gouillart,GaëlVaroquaux

本节介绍使用核心科学模块NumPy和SciPy的基本图像处理和处理。本教程涵盖的一些操作可能对于其他类型的多维数组处理比图像处理有用。特别地,子模块scipy.ndimage提供对n维NumPy数组操作的函数。

也可以看看

有关更高级的图像处理和图像特定程序,请参阅专用于skimage模块的教程Scikit-image: image processing

图像= 2-D数值阵列

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

这里,image == Numpy数组 np.array

本教程中使用的工具

  • numpy:基本数组操作

  • scipyscipy.ndimage子模块专用于图像处理(n维图像)。请参阅文档

    >>> from scipy import ndimage
    

图像处理中的常见任务

  • 输入/输出,显示图像
  • 基本操作:裁剪,翻转,旋转,...
  • 图像过滤:去噪,锐化
  • 图像分割:标记对应于不同对象的像素
  • 分类
  • 特征提取
  • 注册
  • ...

2.6.1. 打开和写入图像文件

将数组写入文件:

from scipy import misc
f = misc.face()
misc.imsave('face.png', f) # uses the Image module (PIL)
import matplotlib.pyplot as plt
plt.imshow(f)
plt.show()
../../_images/face.png

从映像文件创建numpy数组:

>>> from scipy import misc
>>> face = misc.face()
>>> misc.imsave('face.png', face) # First we need to create the PNG file
>>> face = misc.imread('face.png')
>>> type(face)
<... 'numpy.ndarray'>
>>> face.shape, face.dtype
((768, 1024, 3), dtype('uint8'))

对于8位图像,dtype为uint8(0-255)

打开原始文件(摄像机,3-D图像)

>>> face.tofile('face.raw') # Create raw file
>>> face_from_raw = np.fromfile('face.raw', dtype=np.uint8)
>>> face_from_raw.shape
(2359296,)
>>> face_from_raw.shape = (768, 1024, 3)

需要知道图像的形状和类型(如何分离数据字节)。

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

>>> face_memmap = np.memmap('face.raw', dtype=np.uint8, shape=(768, 1024, 3))

(数据从文件读取,而不是加载到内存中)

处理图像文件列表

>>> for i in range(10):
... im = np.random.randint(0, 256, 10000).reshape((100, 100))
... misc.imsave('random_%02d.png' % i, im)
>>> from glob import glob
>>> filelist = glob('random*.png')
>>> filelist.sort()

2.6.2. 显示图像

使用matplotlibimshowmatplotlib

>>> f = misc.face(gray=True)  # retrieve a grayscale image
>>> import matplotlib.pyplot as plt
>>> plt.imshow(f, cmap=plt.cm.gray)
<matplotlib.image.AxesImage object at 0x...>

通过设置最小和最大值增加对比度:

>>> plt.imshow(f, cmap=plt.cm.gray, vmin=30, vmax=200)        
<matplotlib.image.AxesImage object at 0x...>
>>> # Remove axes and ticks
>>> plt.axis('off')
(-0.5, 1023.5, 767.5, -0.5)

绘制轮廓线:

>>> plt.contour(f, [50, 200])        
<matplotlib.contour.QuadContourSet ...>
../../_images/sphx_glr_plot_display_face_001.png

[Python源代码]

要精细检查强度变化,请使用interpolation='nearest'

>>> plt.imshow(f[320:340, 510:530], cmap=plt.cm.gray)        
<matplotlib.image.AxesImage object at 0x...>
>>> plt.imshow(f[320:340, 510:530], cmap=plt.cm.gray, interpolation='nearest')
<matplotlib.image.AxesImage object at 0x...>
../../_images/sphx_glr_plot_interpolation_face_001.png

[Python源代码]

也可以看看

3-D可视化:Mayavi

请参阅3D plotting with Mayavi

  • 图像平面小部件
  • 等面
  • ...
../../_images/decorations.png

2.6.3. 基本操作

图像是数组:使用整个numpy机制。

../../_images/axis_convention.png
>>> face = misc.face(gray=True)
>>> face[0, 40]
127
>>> # Slicing
>>> face[10:13, 20:23]
array([[141, 153, 145],
[133, 134, 125],
[ 96, 92, 94]], dtype=uint8)
>>> face[100:120] = 255
>>>
>>> lx, ly = face.shape
>>> X, Y = np.ogrid[0:lx, 0:ly]
>>> mask = (X - lx / 2) ** 2 + (Y - ly / 2) ** 2 > lx * ly / 4
>>> # Masks
>>> face[mask] = 0
>>> # Fancy indexing
>>> face[range(400), range(400)] = 255
../../_images/sphx_glr_plot_numpy_array_001.png

[Python源代码]

2.6.3.1. 统计信息

>>> face = misc.face(gray=True)
>>> face.mean()
113.48026784261067
>>> face.max(), face.min()
(250, 0)

np.histogram

练习

  • 作为数组打开scikit-image标志(http://scikit-image.org/_static/img/logo.png),或者您拥有的图片你的电脑。
  • 裁剪图像的有意义的部分,例如徽标中的python圆。
  • 使用matplotlib显示图像数组。更改插值方法和缩放以查看差异。
  • 将图像转换为灰度
  • 通过更改图像的最小和最大值来增加图像的对比度。可选:使用scipy.stats.scoreatpercentile(阅读docstring!)以使5%的最暗像素和5%最亮像素饱和。
  • 将数组保存为两种不同的文件格式(png,jpg,tiff)
../../_images/scikit_image_logo.png

2.6.3.2. 几何变换

>>> face = misc.face(gray=True)
>>> lx, ly = face.shape
>>> # Cropping
>>> crop_face = face[lx / 4: - lx / 4, ly / 4: - ly / 4]
>>> # up <-> down flip
>>> flip_ud_face = np.flipud(face)
>>> # rotation
>>> rotate_face = ndimage.rotate(face, 45)
>>> rotate_face_noreshape = ndimage.rotate(face, 45, reshape=False)
../../_images/sphx_glr_plot_geom_face_001.png

[Python源代码]

2.6.4. 图片过滤

局部过滤器:用相邻像素值的函数替换像素值。

邻域:正方形(选择大小),磁盘或更复杂的结构元素

../../_images/kernels.png

2.6.4.1. 模糊/平滑

高斯滤波器来自scipy.ndimage

>>> from scipy import misc
>>> face = misc.face(gray=True)
>>> blurred_face = ndimage.gaussian_filter(face, sigma=3)
>>> very_blurred = ndimage.gaussian_filter(face, sigma=5)

均匀过滤器

>>> local_mean = ndimage.uniform_filter(face, size=11)
../../_images/sphx_glr_plot_blur_001.png

[Python源代码]

2.6.4.2. 锐化

锐化模糊图片:

>>> from scipy import misc
>>> face = misc.face(gray=True).astype(float)
>>> blurred_f = ndimage.gaussian_filter(face, 3)

通过添加拉普拉斯算子的近似来增加边的权重:

>>> filter_blurred_f = ndimage.gaussian_filter(blurred_f, 1)
>>> alpha = 30
>>> sharpened = blurred_f + alpha * (blurred_f - filter_blurred_f)
../../_images/sphx_glr_plot_sharpen_001.png

[Python源代码]

2.6.4.3. 去噪

嘈杂的脸:

>>> from scipy import misc
>>> f = misc.face(gray=True)
>>> f = f[230:290, 220:320]
>>> noisy = f + 0.4 * f.std() * np.random.random(f.shape)

高斯滤波器也平滑噪声输出...和边缘:

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

大多数局部线性各向同性滤波器模糊图像(ndimage.uniform_filter

中值滤波器可更好地保留边缘:

>>> med_denoised = ndimage.median_filter(noisy, 3)
../../_images/sphx_glr_plot_face_denoise_001.png

[Python源代码]

中值滤波器:对于直边界(低曲率)的结果更好:

>>> im = np.zeros((20, 20))
>>> im[5:-5, 5:-5] = 1
>>> im = ndimage.distance_transform_bf(im)
>>> im_noise = im + 0.2 * np.random.randn(*im.shape)
>>> im_med = ndimage.median_filter(im_noise, 3)
../../_images/sphx_glr_plot_denoising_001.png

[Python源代码]

其他排序过滤器:ndimage.maximum_filterndimage.percentile_filter

其他本地非线性过滤器:Wiener(scipy.signal.wiener)等。

非本地过滤器

练习:去噪

  • 创建具有多个对象(圆,椭圆,正方形或随机形状)的二进制图像(0和1)。
  • 添加一些噪声(例如,20%的噪声)
  • 尝试两种不同的去噪方法去噪图像:高斯滤波和中值滤波。
  • 比较两个不同的去噪图像的直方图。哪一个是最接近原始(无噪声)图像的直方图?

也可以看看

更多去噪滤波器可在skimage.denoising中获得,请参阅Scikit-image: image processing教程。

2.6.4.4. 数学形态

有关数学形态学的定义,请参见维基百科

探测具有简单形状的图像(结构元素),并根据形状在本地适合或错过图像来修改此图像。

结构元素

>>> el = ndimage.generate_binary_structure(2, 1)
>>> el
array([[False, True, False],
[ True, True, True],
[False, True, False]], dtype=bool)
>>> el.astype(np.int)
array([[0, 1, 0],
[1, 1, 1],
[0, 1, 0]])
../../_images/diamond_kernel.png

侵蚀 =最小过滤器。将像素的值替换为结构化元素覆盖的最小值。

>>> a = np.zeros((7,7), dtype=np.int)
>>> a[1:6, 2:5] = 1
>>> a
array([[0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 1, 1, 0, 0],
[0, 0, 1, 1, 1, 0, 0],
[0, 0, 1, 1, 1, 0, 0],
[0, 0, 1, 1, 1, 0, 0],
[0, 0, 1, 1, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0]])
>>> ndimage.binary_erosion(a).astype(a.dtype)
array([[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0]])
>>> #Erosion removes objects smaller than the structure
>>> ndimage.binary_erosion(a, structure=np.ones((5,5))).astype(a.dtype)
array([[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0]])
../../_images/morpho_mat.png

扩张:最大过滤器:

>>> a = np.zeros((5, 5))
>>> a[2, 2] = 1
>>> a
array([[ 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0.],
[ 0., 0., 1., 0., 0.],
[ 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0.]])
>>> ndimage.binary_dilation(a).astype(a.dtype)
array([[ 0., 0., 0., 0., 0.],
[ 0., 0., 1., 0., 0.],
[ 0., 1., 1., 1., 0.],
[ 0., 0., 1., 0., 0.],
[ 0., 0., 0., 0., 0.]])

也适用于灰值图像:

>>> np.random.seed(2)
>>> im = np.zeros((64, 64))
>>> x, y = (63*np.random.random((2, 8))).astype(np.int)
>>> im[x, y] = np.arange(8)
>>> bigger_points = ndimage.grey_dilation(im, size=(5, 5), structure=np.ones((5, 5)))
>>> square = np.zeros((16, 16))
>>> square[4:-4, 4:-4] = 1
>>> dist = ndimage.distance_transform_bf(square)
>>> dilate_dist = ndimage.grey_dilation(dist, size=(3, 3), \
... structure=np.ones((3, 3)))
../../_images/sphx_glr_plot_greyscale_dilation_001.png

[Python源代码]

打开:erosion + dilation:

>>> a = np.zeros((5,5), dtype=np.int)
>>> a[1:4, 1:4] = 1; a[4, 4] = 1
>>> a
array([[0, 0, 0, 0, 0],
[0, 1, 1, 1, 0],
[0, 1, 1, 1, 0],
[0, 1, 1, 1, 0],
[0, 0, 0, 0, 1]])
>>> # Opening removes small objects
>>> ndimage.binary_opening(a, structure=np.ones((3,3))).astype(np.int)
array([[0, 0, 0, 0, 0],
[0, 1, 1, 1, 0],
[0, 1, 1, 1, 0],
[0, 1, 1, 1, 0],
[0, 0, 0, 0, 0]])
>>> # Opening can also smooth corners
>>> ndimage.binary_opening(a).astype(np.int)
array([[0, 0, 0, 0, 0],
[0, 0, 1, 0, 0],
[0, 1, 1, 1, 0],
[0, 0, 1, 0, 0],
[0, 0, 0, 0, 0]])

应用:remove noise:

>>> square = np.zeros((32, 32))
>>> square[10:-10, 10:-10] = 1
>>> np.random.seed(2)
>>> x, y = (32*np.random.random((2, 20))).astype(np.int)
>>> square[x, y] = 1
>>> open_square = ndimage.binary_opening(square)
>>> eroded_square = ndimage.binary_erosion(square)
>>> reconstruction = ndimage.binary_propagation(eroded_square, mask=square)
../../_images/sphx_glr_plot_propagation_001.png

[Python源代码]

关闭:扩张+侵蚀

许多其他数学形态运算:命中和遗漏变换,tophat等。

2.6.5. 特征提取

2.6.5.1. 边缘检测

合成数据:

>>> im = np.zeros((256, 256))
>>> im[64:-64, 64:-64] = 1
>>>
>>> im = ndimage.rotate(im, 15, mode='constant')
>>> im = ndimage.gaussian_filter(im, 8)

使用渐变运算符Sobel)查找高强度变化:

>>> sx = ndimage.sobel(im, axis=0, mode='constant')
>>> sy = ndimage.sobel(im, axis=1, mode='constant')
>>> sob = np.hypot(sx, sy)
../../_images/sphx_glr_plot_find_edges_001.png

[Python源代码]

2.6.5.2. 分段

  • 基于直方图的分割(无空间信息)
>>> n = 10
>>> l = 256
>>> im = np.zeros((l, l))
>>> np.random.seed(1)
>>> points = l*np.random.random((2, n**2))
>>> im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
>>> im = ndimage.gaussian_filter(im, sigma=l/(4.*n))
>>> mask = (im > im.mean()).astype(np.float)
>>> mask += 0.1 * im
>>> img = mask + 0.2*np.random.randn(*mask.shape)
>>> hist, bin_edges = np.histogram(img, bins=60)
>>> bin_centers = 0.5*(bin_edges[:-1] + bin_edges[1:])
>>> binary_img = img > 0.5
../../_images/sphx_glr_plot_histo_segmentation_001.png

[Python源代码]

使用数学形态学来清除结果:

>>> # Remove small white regions
>>> open_img = ndimage.binary_opening(binary_img)
>>> # Remove small black hole
>>> close_img = ndimage.binary_closing(open_img)
../../_images/sphx_glr_plot_clean_morpho_001.png

[Python源代码]

练习

检查重建操作(侵蚀+传播)是否产生比打开/关闭更好的结果:

>>> eroded_img = ndimage.binary_erosion(binary_img)
>>> reconstruct_img = ndimage.binary_propagation(eroded_img, mask=binary_img)
>>> tmp = np.logical_not(reconstruct_img)
>>> eroded_tmp = ndimage.binary_erosion(tmp)
>>> reconstruct_final = np.logical_not(ndimage.binary_propagation(eroded_tmp, mask=tmp))
>>> np.abs(mask - close_img).mean()
0.00727836...
>>> np.abs(mask - reconstruct_final).mean()
0.00059502...

练习

检查第一去噪步骤(例如,使用中值滤波器)如何修改直方图,并检查所得的基于直方图的分割是否更准确。

也可以看看

更高级的分割算法在scikit-image中找到:见Scikit-image: image processing

也可以看看

其他科学包提供了可用于图像处理的算法。在这个例子中,我们使用scikit-learn的光谱聚类函数来分割粘贴的对象。

>>> from sklearn.feature_extraction import image
>>> from sklearn.cluster import spectral_clustering
>>> l = 100
>>> x, y = np.indices((l, l))
>>> center1 = (28, 24)
>>> center2 = (40, 50)
>>> center3 = (67, 58)
>>> center4 = (24, 70)
>>> radius1, radius2, radius3, radius4 = 16, 14, 15, 14
>>> circle1 = (x - center1[0])**2 + (y - center1[1])**2 < radius1**2
>>> circle2 = (x - center2[0])**2 + (y - center2[1])**2 < radius2**2
>>> circle3 = (x - center3[0])**2 + (y - center3[1])**2 < radius3**2
>>> circle4 = (x - center4[0])**2 + (y - center4[1])**2 < radius4**2
>>> # 4 circles
>>> img = circle1 + circle2 + circle3 + circle4
>>> mask = img.astype(bool)
>>> img = img.astype(float)
>>> img += 1 + 0.2*np.random.randn(*img.shape)
>>> # Convert the image into a graph with the value of the gradient on
>>> # the edges.
>>> graph = image.img_to_graph(img, mask=mask)
>>> # Take a decreasing function of the gradient: we take it weakly
>>> # dependant from the gradient the segmentation is close to a voronoi
>>> graph.data = np.exp(-graph.data/graph.data.std())
>>> labels = spectral_clustering(graph, n_clusters=4, eigen_solver='arpack')
>>> label_im = -np.ones(mask.shape)
>>> label_im[mask] = labels
../../_images/image_spectral_clustering.png

2.6.6. 测量对象属性:ndimage.measurements

合成数据:

>>> n = 10
>>> l = 256
>>> im = np.zeros((l, l))
>>> points = l*np.random.random((2, n**2))
>>> im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
>>> im = ndimage.gaussian_filter(im, sigma=l/(4.*n))
>>> mask = im > im.mean()
  • 连接组件的分析

标签连接的组件:ndimage.label

>>> label_im, nb_labels = ndimage.label(mask)
>>> nb_labels # how many regions?
16
>>> plt.imshow(label_im)
<matplotlib.image.AxesImage object at 0x...>
../../_images/sphx_glr_plot_synthetic_data_001.png

[Python源代码]

计算大小,平均值等。的区域:

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

清理小连接组件:

>>> mask_size = sizes < 1000
>>> remove_pixel = mask_size[label_im]
>>> remove_pixel.shape
(256, 256)
>>> label_im[remove_pixel] = 0
>>> plt.imshow(label_im)
<matplotlib.image.AxesImage object at 0x...>

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

>>> labels = np.unique(label_im)
>>> label_im = np.searchsorted(labels, label_im)
../../_images/sphx_glr_plot_measure_data_001.png

[Python源代码]

查找包围物体的感兴趣区域:

>>> slice_x, slice_y = ndimage.find_objects(label_im==4)[0]
>>> roi = im[slice_x, slice_y]
>>> plt.imshow(roi)
<matplotlib.image.AxesImage object at 0x...>
../../_images/sphx_glr_plot_find_object_001.png

[Python源代码]

其他空间尺寸:ndimage.center_of_massndimage.maximum_position等。

可以在有限的分割应用范围外使用。

示例:block mean:

>>> from scipy import misc
>>> f = misc.face(gray=True)
>>> sx, sy = f.shape
>>> X, Y = np.ogrid[0:sx, 0:sy]
>>> regions = (sy//6) * (X//4) + (Y//6) # note that we use broadcasting
>>> block_mean = ndimage.mean(f, labels=regions, index=np.arange(1,
... regions.max() +1))
>>> block_mean.shape = (sx // 4, sy // 6)
../../_images/sphx_glr_plot_block_mean_001.png

[Python源代码]

当区域是常规块时,使用步幅技巧(Example: fake dimensions with strides)更为有效。

非规则间隔块:径向平均值:

>>> sx, sy = f.shape
>>> X, Y = np.ogrid[0:sx, 0:sy]
>>> r = np.hypot(X - sx/2, Y - sy/2)
>>> rbin = (20* r/r.max()).astype(np.int)
>>> radial_mean = ndimage.mean(f, labels=rbin, index=np.arange(1, rbin.max() +1))
../../_images/sphx_glr_plot_radial_mean_001.png

[Python源代码]

  • 其他措施

相关函数,傅立叶/小波谱等

数学形态学的一个例子:粒度

>>> def disk_structure(n):
... struct = np.zeros((2 * n + 1, 2 * n + 1))
... x, y = np.indices((2 * n + 1, 2 * n + 1))
... mask = (x - n)**2 + (y - n)**2 <= n**2
... struct[mask] = 1
... return struct.astype(np.bool)
...
>>>
>>> def granulometry(data, sizes=None):
... s = max(data.shape)
... if sizes == None:
... sizes = range(1, s/2, 2)
... granulo = [ndimage.binary_opening(data, \
... structure=disk_structure(n)).sum() for n in sizes]
... return granulo
...
>>>
>>> np.random.seed(1)
>>> n = 10
>>> l = 256
>>> im = np.zeros((l, l))
>>> points = l*np.random.random((2, n**2))
>>> im[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
>>> im = ndimage.gaussian_filter(im, sigma=l/(4.*n))
>>>
>>> mask = im > im.mean()
>>>
>>> granulo = granulometry(mask, sizes=np.arange(2, 19, 4))
../../_images/sphx_glr_plot_granulo_001.png

[Python源代码]

也可以看看

更多关于图像处理: