这个教程中使用的工具:
from scipy import ndimage
- 一些例子用到了使用np.array的特殊的工具箱:
图像中的常见问题有:
- 输入/输出:呈现图像
- 基本操作:裁剪、翻转、旋转、比例缩放……
- 图像滤镜:消噪,锐化
- 图像分割:不同对应对象的像素标记
更有力和完整的模块:
- OpenCV (Python绑定)
- CellProfiler
- ITK,Python绑定
- 更多……
1、打开和读写图像文件
将一个数组写入文件:
In [1]: from scipy import misc
In [2]: l = misc.lena()
In [3]: misc.imsave('lena.png', l) # uses the Image module (PIL)
In [4]: import pylab as pl
In [5]: pl.imshow(l)
Out[5]: <matplotlib.image.AxesImage at 0x4118110>
从一个图像文件创建数组:
import matplotlib.pyplot as plt
from scipy import misc
import pylab as pl
img1=misc.imread('Lenna.png',mode='L')
img2=misc.imread('Lenna.png',mode='RGB')
img3=misc.imread('Lenna.png',mode='RGBA')
pl.imshow(img2,cmap='gray')
pl.show()
misc.imsave('lena.png',img2)
或
from scipy import ndimage
ndimage.imread('lena.png',mode='L')
ndimage.imread('lena.png', mode='RGB')
ndimage.imread('lena.png', mode='RGBA')
打开一个raw文件(相机, 3-D图像)
from scipy import misc
import numpy as np
import pylab as pl
img=misc.imread('Lenna.png',mode='L')
img.tofile('lena.raw')
lena_from_raw = np.fromfile('lena.raw',np.uint8)
print(lena_from_raw.shape)
lena_from_raw.shape = (512, 512)
pl.imshow(lena_from_raw,cmap='gray')
pl.show()
import os
os.remove('lena.raw')
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
对于大数据,使用np.memmap进行内存映射:
from scipy import misc
import numpy as np
import pylab as pl
img=misc.imread('Lenna.png',mode='RGB')
img.tofile('lena.raw')
lena_memmap = np.memmap('lena.raw', dtype=np.uint8, shape=(512,512,3))
pl.imshow(lena_memmap,cmap='gray')
pl.show()
import os
os.remove('lena.raw')
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
批量处理图像文件:
In [22]: for i in range(10):
....: im = np.random.random_integers(0, 255, 10000).reshape((100, 100))
....: misc.imsave(
....:
In [23]: from glob import glob
In [24]: filelist = glob(
In [25]: filelist.sort()
2、呈现图像
from scipy import misc
import numpy as np
import pylab as pl
import matplotlib.pyplot as plt
img=misc.imread('Lenna.png',mode='L')
plt.imshow(img,cmap=plt.cm.gray)
plt.imshow(img, cmap='gray', vmin=30, vmax=200)
plt.imshow(img[200:220, 200:220], cmap=plt.cm.gray, interpolation='nearest')
plt.axis('off')
plt.contour(img, [60, 211])
plt.show()
pl.imshow(img,cmap=plt.cm.gray)
pl.imshow(img, cmap='gray', vmin=30, vmax=200)
pl.imshow(img[200:220, 200:220], cmap=pl.cm.gray, interpolation='nearest')
pl.axis('off')
pl.contour(img, [60, 211])
pl.show()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
其它包有时使用图形工具箱来可视化(GTK,Qt)
from scipy import misc
import numpy as np
import skimage.io as im_io
img=misc.imread('Lenna.png',mode='L')
im_io.use_plugin('gtk', 'imshow')
im_io.imshow(img)
im_io.show()
3-D可视化:Mayavi
参见可用Mayavi进行3-D绘图和体积数据
3、基本操作
图像是数组(numpy array)。
>>> lena = misc.lena()
>>> lena[0, 40]
166
>>>
>>> lena[10:13, 20:23]
array([[158, 156, 157],
[157, 155, 155],
[157, 157, 158]])
>>> lena[100:120] = 255
>>>
>>> lx, ly = lena.shape
>>> X, Y = np.ogrid[0:lx, 0:ly]
>>> mask = (X - lx/2)**2 + (Y - ly/2)**2 > lx*ly/4
>>>
>>> lena[mask] = 0
>>>
>>> lena[range(400), range(400)] = 255
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
统计信息
from scipy import misc
import numpy as np
import pylab as pl
import matplotlib.pyplot as plt
import cv2
lena=misc.imread('Lenna.png',mode='L')
lena.mean()
lena.max()
lena.min()
hist, bins = np.histogram(lena.ravel(), bins=256, range=[0,256],normed=True)
plt.plot(hist)
plt.legend(loc='best')
plt.show()
plt.hist(lena.ravel(),256,[0,256]); plt.show()
hist_full = cv2.calcHist([lena],[0],None,[256],[0,256])
plt.plot(hist_full)
plt.xlim([0,256])
plt.show()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
几何转换
from scipy import misc
import numpy as np
import scipy
import pylab as pl
import matplotlib.pyplot as plt
import cv2
lena=misc.imread('Lenna.png',mode='L')
lx, ly = lena.shape
crop_lena = lena[lx/4:-lx/4, ly/4:-ly/4]
flip_ud_lena = np.flipud(lena)
flip_lr_lena = np.fliplr(lena)
rotate_lena = ndimage.rotate(lena, 45)
rotate_lena_noreshape = ndimage.rotate(lena, 45, reshape=False)
zoomed_lena = ndimage.zoom(lena, 2)
zoomed_lena2 = ndimage.zoom(lena, 0.5)
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
4、图像滤镜
局部滤镜
:用相邻像素值的函数替代当前像素的值。
相邻:方形(指定大小),圆形, 或者更多复杂的结构元素。
模糊/平滑
scipy.ndimage中的高斯滤镜:
>>> from scipy import misc
>>> from scipy import ndimage
>>> lena = misc.lena()
>>> blurred_lena = ndimage.gaussian_filter(lena, sigma=3)
>>> very_blurred = ndimage.gaussian_filter(lena, sigma=5)
均匀滤镜
>>> local_mean = ndimage.uniform_filter(lena, size=11)
锐化
锐化模糊图像:
from scipy import misc
import pylab as pl
from scipy import ndimage
lena=misc.imread('Lenna.png',mode='L')
blurred_l = ndimage.gaussian_filter(lena, 3)
filter_blurred_l = ndimage.gaussian_filter(blurred_l, 1)
alpha = 30
sharpened = blurred_l + alpha * (blurred_l - filter_blurred_l)
pl.subplot(131),pl.imshow(blurred_l,cmap='gray'),pl.axis('off'),pl.title('blurred_l')
pl.subplot(132),pl.imshow(filter_blurred_l,cmap='gray'),pl.axis('off'),pl.title('filter_blurred_l')
pl.subplot(133),pl.imshow(sharpened,cmap='gray'),pl.axis('off'),pl.title('sharpened')
pl.show()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
消噪
from scipy import misc
import pylab as pl
from scipy import ndimage
import numpy as np
lena=misc.imread('Lenna.png',mode='L')
l = lena[230:310, 210:350]
noisy = l + 0.4*l.std()*np.random.random(l.shape)
gauss_denoised = ndimage.gaussian_filter(noisy, 2)
unif_denoised=ndimage.uniform_filter(noisy,2)
med_denoised = ndimage.median_filter(noisy, 3)
pl.subplot(141),pl.imshow(noisy,cmap='gray'),pl.axis('off'),pl.title('noisy')
pl.subplot(142),pl.imshow(gauss_denoised,cmap='gray'),pl.axis('off'),pl.title('gauss_denoised')
pl.subplot(143),pl.imshow(unif_denoised,cmap='gray'),pl.axis('off'),pl.title('unif_denoised')
pl.subplot(144),pl.imshow(med_denoised,cmap='gray'),pl.axis('off'),pl.title('med_denoised')
pl.show()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
其它排序滤波器:
ndimage.maximum_filter,ndimage.percentile_filter
其它局部非线性滤波器:
维纳滤波器(scipy.signal.wiener)等
非局部滤波器
总变差(TV)消噪。找到新的图像让图像的总变差(正态L1梯度的积分)变得最小,当接近测量图像时
from scipy import misc
import pylab as pl
from scipy import ndimage
import numpy as np
from skimage.restoration import denoise_tv_bregman,denoise_tv_chambolle
lena=misc.imread('Lenna.png',mode='L')
l = lena[230:310, 210:350]
noisy = l + 0.4*l.std()*np.random.random(l.shape)
tv_denoised_b = denoise_tv_bregman(noisy, weight=10)
tv_denoised_c = denoise_tv_chambolle(noisy, weight=10)
ax=pl.figure()
ax.add_subplot(131),pl.imshow(noisy,cmap='gray'),pl.axis('off'),pl.title('noisy')
pl.show()
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
5、数学形态学
结构元素
>>> 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]])
腐蚀 = 最小化滤镜。用结构元素覆盖的像素的最小值替代一个像素值:
>>> 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]])
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
膨胀:最大化滤镜:
>>> 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)
>>> 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)))
开操作:腐蚀+膨胀
应用:移除噪声
>>> 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)
示例代码
闭操作:膨胀+腐蚀
许多其它数学分形:击中(hit)和击不中(miss)变换,tophat等等。
6、特征提取
边缘检测
使用梯度操作(Sobel)来找到搞强度的变化
示例源码
canny滤镜
Canny滤镜可以从skimage中获取(文档)
>>>
>>>
>>> im += 0.1*np.random.random(im.shape)
>>> edges = canny(im, 1, 0.4, 0.2)
>>> edges = canny(im, 3, 0.3, 0.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
示例代码
自动阈值:使用高斯混合模型
>>> mask = (im > im.mean()).astype(np.float)
>>> mask += 0.1 * im
>>> img = mask + 0.3*np.random.randn(*mask.shape)
>>> from sklearn.mixture import GMM
>>> classif = GMM(n_components=2)
>>> classif.fit(img.reshape((img.size, 1)))
GMM(...)
>>> classif.means_
array([[ 0.9353155 ],
[-0.02966039]])
>>> np.sqrt(classif.covars_).ravel()
array([ 0.35074631, 0.28225327])
>>> classif.weights_
array([ 0.40989799, 0.59010201])
>>> threshold = np.mean(classif.means_)
>>> binary_img = img > threshold
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
使用数学形态学来清理结果:
>>>
>>> open_img = ndimage.binary_opening(binary_img)
>>>
>>> close_img = ndimage.binary_closing(open_img)