In [2]:
%matplotlib inline
import numpy as np
作者：Emmanuelle Gouillart, Gaël Varoquaux
原文: http://www.scipy-lectures.org/advanced/image_processing/index.html
这个部分解决用核心的科学模块NumPy和SciPy做基本的图像操作和处理。这个教程中涵盖的一些操作可能对于一些其他类型的多维度数据处理比对图像处理更加有用。特别是，子摸块scipy.ndimage提供了在N维Numpy数组上操作的方法。
也看一下: 对于更高级的图像处理和图像特有的程序，见专注于skimage模块教程Scikit-image: 图像处理。
图像 = 2-D 数值数组
(或者 3-D: CT, MRI, 2D + time; 4-D, ...)
这里, 图像 == Numpy 数组 np.array
本教程中使用的工具:
numpy: 基础的数组操作scipy:scipy.ndimage专注于图像处理的子模块 (n维 图像)。见文档:
In [3]:
from scipy import ndimage
图像处理中的常见任务:
- 输入/输出、显示图像
- 基础操作: 剪切, 翻转、旋转...
- 图像过滤: 降噪, 锐化
- 图形分割: 根据不同的对象标记像素
- 分类
- 特征提取
- 配准
- ...
章节内容
打开和写入图像文件
显示图像
基础操作
统计信息
几何图像变换
图像过滤
模糊/光滑
锐化
降噪
数学形态学
特征提取
边缘检测
分隔
测量对象属性: ndimage.measurements
2.6.1 打开和写入图像文件¶
将数组写入文件:
In [4]:
from scipy import misc
f = misc.face()
misc.imsave('face.png', f) # 使用图像模块 (PIL)
import matplotlib.pyplot as plt
plt.imshow(f)
plt.show()
从图像文件创建一个numpy数组:
In [5]:
from scipy import misc
face = misc.face()
misc.imsave('face.png', face) # 首先我们需要创建这个PNG文件
face = misc.imread('face.png')
type(face)
Out[5]:
In [6]:
face.shape, face.dtype
Out[6]:
对于8位的图像 (0-255) dtype是uint8
打开raw文件 (照相机, 3-D 图像)
In [7]:
face.tofile('face.raw') # 创建raw文件
face_from_raw = np.fromfile('face.raw', dtype=np.uint8)
face_from_raw.shape
Out[7]:
In [8]:
face_from_raw.shape = (768, 1024, 3)
需要知道图像的shape和dtype (如何去分离数据类型)。
对于大数据, 使用np.memmap来做内存映射:
In [9]:
face_memmap = np.memmap('face.raw', dtype=np.uint8, shape=(768, 1024, 3))
(数据从文件中读取，并没有加载到内存)
处理一组图像文件
In [10]:
for i in range(10):
im = np.random.random_integers(0, 255, 10000).reshape((100, 100))
misc.imsave('random_%02d.png' % i, im)
from glob import glob
filelist = glob('random*.png')
filelist.sort()
2.6.2 显示图像¶
使用matplotlib和imshow在matplotlib图形内部显示图像:
In [11]:
f = misc.face(gray=True) # 取回灰度图像
import matplotlib.pyplot as plt
plt.imshow(f, cmap=plt.cm.gray)
Out[11]:
通过设置最小和最大值增加对比度:
In [12]:
plt.imshow(f, cmap=plt.cm.gray, vmin=30, vmax=200)
Out[12]:
In [13]:
plt.imshow(f, cmap=plt.cm.gray, vmin=30, vmax=200)
# 删除座标轴和刻度
plt.axis('off')
Out[13]:
画出轮廓线:
In [14]:
plt.imshow(f, cmap=plt.cm.gray, vmin=30, vmax=200)
# 删除座标轴和刻度
plt.axis('off')
plt.contour(f, [50, 200])
Out[14]:
对于要精确检查的密度变量，使用interpolation='nearest':
In [15]:
plt.imshow(f[320:340, 510:530], cmap=plt.cm.gray)
Out[15]:
In [16]:
plt.imshow(f[320:340, 510:530], cmap=plt.cm.gray, interpolation='nearest')
Out[16]:
In [17]:
face = misc.face(gray=True)
face[0, 40]
Out[17]:
In [18]:
# 切片
face[10:13, 20:23]
Out[18]:
In [19]:
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
2.6.3.1 统计信息¶
In [20]:
face = misc.face(gray=True)
face.mean()
Out[20]:
In [21]:
face.max(), face.min()
Out[21]:
np.histogram
练习
- 将
scikit-imagelogo作为数组打开 (http://scikit-image.org/_static/img/logo.png), 或者在你电脑上的其他图像。 - 剪切图像的有意义部分，例如，logo中的python圆形。
- 用
matplotlib显示图像数组。改变interpolation方法并且放大看一下差异。 - 将你的图像改变为灰度
- 通过改变它的最小最大值增加图像的对比度。选做：使用
scipy.stats.scoreatpercentile(读一下文本字符串!) 来饱和最黑5%的像素和最亮5%的像素。 - 将数组保存为两个不同的文件格式 (png, jpg, tiff)
2.6.3.2 几何图像变换¶
In [22]:
face = misc.face(gray=True)
lx, ly = face.shape
# 剪切
crop_face = face[lx / 4: - lx / 4, ly / 4: - ly / 4]
# up <-> down 翻转
flip_ud_face = np.flipud(face)
# 旋转
rotate_face = ndimage.rotate(face, 45)
rotate_face_noreshape = ndimage.rotate(face, 45, reshape=False)
2.6.4.1 模糊 / 光滑¶
高斯过滤器 来自 scipy.ndimage:
In [23]:
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)
均匀过滤器
In [24]:
local_mean = ndimage.uniform_filter(face, size=11)
In [25]:
from scipy import misc
face = misc.face(gray=True).astype(float)
blurred_f = ndimage.gaussian_filter(face, 3)
通过添加Laplacian的近似值来增加边缘的权重:
In [26]:
filter_blurred_f = ndimage.gaussian_filter(blurred_f, 1)
alpha = 30
sharpened = blurred_f + alpha * (blurred_f - filter_blurred_f)
2.6.4.3 降噪¶
有噪音的脸:
In [27]:
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)
高斯过滤器光滑了噪音... 以及边缘:
In [28]:
gauss_denoised = ndimage.gaussian_filter(noisy, 2)
绝大多数局部线性各向同性过滤器模糊图像 (ndimage.uniform_filter)
中位数过滤器更好的保留的边缘:
In [29]:
med_denoised = ndimage.median_filter(noisy, 3)
中位数过滤器: 对直边缘的结果更好 (低曲度):
In [30]:
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)
其他排名过滤器：ndimage.maximum_filter、ndimage.percentile_filter
其他的局部非线性过滤器: Wiener (scipy.signal.wiener) 等等.
非局部过滤器
练习: 降噪
- 创建一个带有一些对象 (圆形、椭圆、正方形或随机形状) 的二元图像 (0 和 1) .
- 添加一些噪音 (例如, 20% 的噪音)
- 用两种不同的降噪方法来降噪这个图像: 高斯过滤器和中位数过滤器。
- 比较两种不同降噪方法的直方图。哪一个与原始图像（无噪音）的直方图最接近?
也看一下： 在skimage.denoising中有更多的的降噪过滤器可用，见教程Scikit-image: 图像处理.
2.6.4.4 数学形态学¶
看一下wikipedia上的数学形态学定义。
探索一个简单形状 (结构化的元素) 的图像，然后根据这个形状是如何局部适应或不适应这个图像来修改这个图像。
结构化元素:
In [31]:
el = ndimage.generate_binary_structure(2, 1)
el
Out[31]:
In [32]:
el.astype(np.int)
Out[32]:
风化(Erosion) = 最小值过滤器。用结构化元素所覆盖的最小值来替换像素的值。:
In [33]:
a = np.zeros((7,7), dtype=np.int)
a[1:6, 2:5] = 1
a
Out[33]:
In [34]:
ndimage.binary_erosion(a).astype(a.dtype)
Out[34]:
In [35]:
#风化删除了比结构小的对象
ndimage.binary_erosion(a, structure=np.ones((5,5))).astype(a.dtype)
Out[35]:
扩大(Dilation)：最大值过滤器：
In [36]:
a = np.zeros((5, 5))
a[2, 2] = 1
a
Out[36]:
In [37]:
ndimage.binary_dilation(a).astype(a.dtype)
Out[37]:
对于灰度值图像同样适用:
In [38]:
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)))
Opening: erosion + dilation:
In [39]:
a = np.zeros((5,5), dtype=np.int)
a[1:4, 1:4] = 1; a[4, 4] = 1
a
Out[39]:
In [40]:
# Opening 删除小对象
ndimage.binary_opening(a, structure=np.ones((3,3))).astype(np.int)
Out[40]:
In [41]:
# Opening 也可以光滑转角
ndimage.binary_opening(a).astype(np.int)
Out[41]:
应用: 去除噪音:
In [42]:
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)
2.6.5.2 分割¶
- 直方图分割 (没有空间信息)
In [43]:
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