使用Numpy和Scipy處理圖像

jopen 10年前發布 | 118K 次閱讀圖形/圖像處理 NumPy
Image manipulation and processing using Numpy and Scipy

翻譯自：http://scipy-lectures.github.com/advanced/image_processing/index.html
作者：Emmanuelle Gouillart, Ga?l Varoquaux
圖像 = 2-D 數值數組
(或者 3-D: CT, MRI, 2D + 時間; 4-D, ...)
這里 圖像 == Numpy數組 np.array</pre> 
這個教程中使用的工具：



numpy：基本數組操作 
</li>

scipy：scipy.ndimage子模塊致力于圖像處理(n維圖像)。參見http://docs.scipy.org/doc/scipy/reference/tutorial/ndimage.html 

from scipy import ndimage
</li>

一些例子用到了使用np.array的特殊的工具箱： 
</li>


Scikit Image 
</li>

scikit-learn 
</li>
</ul>
</ul>
圖像中的常見問題有：



輸入/輸出，呈現圖像 
</li>

基本操作：裁剪、翻轉、旋轉…… 
</li>

圖像濾鏡：消噪，銳化 
</li>

圖像分割：不同對應對象的像素標記 
</li>
</ul>
更有力和完整的模塊：



OpenCV (Python綁定) 
</li>

CellProfiler 
</li>

ITK,Python綁定 
</li>

更多…… 
</li>
</ul>
目錄 



打開和讀寫圖像文件 
</li>

呈現圖像 
</li>

基本操作 
</li>


統計信息 
</li>

幾何轉換 
</li>
</ul>

圖像濾鏡 
</li>


模糊/平滑 
</li>

銳化 
</li>

消噪 
</li>

數學形態學 
</li>
</ul>

特征提取 
</li>


邊緣檢測 
</li>

分割 
</li>
</ul>

測量對象屬性：ndimage.measurements 
</li>

Footnotes 
</li>
</ul>

打開和讀寫圖像文件


將一個數組寫入文件：

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></pre> 
從一個圖像文件創建數組：

In [7]: lena = misc.imread('lena.png')
In [8]: type(lena)
Out[8]: numpy.ndarray
In [9]: lena.shape, lena.dtype
Out[9]: ((512, 512), dtype('uint8'))</pre> 
8位圖像(0-255)的dtype是uint8

打開一個raw文件(相機， 3-D圖像)

In [10]: l.tofile('lena.raw')  # 創建一個raw文件
In [14]: lena_from_raw = np.fromfile('lena.raw', dtype=np.int64)
In [15]: lena_from_raw.shape
Out[15]: (262144,)
In [16]: lena_from_raw.shape = (512, 512)
In [17]: import os
In [18]: os.remove('lena.raw')</pre> 
需要知道圖像的shape和dtype(如何區分隔數據字節)

對于大數據，使用np.memmap進行內存映射：

In [21]: lena_memmap = np.memmap('lena.raw', dtype=np.int64, shape=(512,512))
 (數據從文件讀取，而不是載入內存)

處理一個列表的圖像文件：

In [22]: for i in range(10):
   ....:     im = np.random.randomintegers(0, 255, 10000).reshape((100, 100))
   ....:     misc.imsave('random%02d.png' % i, im)
   ....:     
In [23]: from glob import glob
In [24]: filelist = glob('random*.png')
In [25]: filelist.sort()</pre> 

呈現圖像


使用matplotlib和imshow將圖像呈現在matplotlib圖像(figure)中：

In [29]: l = misc.lena()
In [30]: import matplotlib.pyplot as plt
In [31]: plt.imshow(l, cmap=plt.cm.gray)
Out[31]: <matplotlib.image.AxesImage at 0x4964990></pre> 
通過設置最大最小之增加對比：

In [33]: plt.imshow(l, cmap=plt.cm.gray, vmin=30, vmax=200)
Out[33]: <matplotlib.image.AxesImage at 0x50cb790>
In [34]: plt.axis('off')  # 移除axes和ticks
Out[34]: (-0.5, 511.5, 511.5, -0.5)</pre> 
繪制等高線：¹ 

ln[7]: plt.contour(l, [60, 211])
 更好地觀察強度變化，使用interpolate=‘nearest’：

In [7]: plt.imshow(l[200:220, 200:220], cmap=plt.cm.gray)
Out[7]: <matplotlib.image.AxesImage at 0x3bbe610>
In [8]: plt.imshow(l[200:220, 200:220], cmap=plt.cm.gray, interpolation='nearest')
Out[8]: <matplotlib.image.AxesImage at 0x3ed3250></pre> 
其它包有時使用圖形工具箱來可視化(GTK，Qt)：² 

In [9]: import skimage.io as im_io
In [21]: im_io.use_plugin('gtk', 'imshow')
In [22]: im_io.imshow(l)</pre> 
3-D可視化：Mayavi 

參見可用Mayavi進行3-D繪圖和體積數據 



圖形平面工具 
</li>

等值面 
</li>

…… 
</li>
</ul>

基本操作


圖像是數組：使用整個numpy機理。

 

>>> lena = misc.lena()
>>> lena[0, 40]
166
>>> # Slicing
>>> 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
>>> # Masks
>>> lena[mask] = 0
>>> # Fancy indexing
>>> lena[range(400), range(400)] = 255
 

統計信息


>>> lena = scipy.lena()
>>> lena.mean()
124.04678344726562
>>> lena.max(), lena.min()
(245, 25)
 np.histogram

幾何轉換


>>> lena = scipy.lena()
>>> lx, ly = lena.shape
>>> # Cropping
>>> crop_lena = lena[lx/4:-lx/4, ly/4:-ly/4]
>>> # up <-> down flip
>>> flip_ud_lena = np.flipud(lena)
>>> # rotation
>>> rotate_lena = ndimage.rotate(lena, 45)
>>> rotate_lena_noreshape = ndimage.rotate(lena, 45, reshape=False)
  
示例源碼 

圖像濾鏡


局部濾鏡：用相鄰像素值的函數替代當前像素的值。
相鄰：方形(指定大小)，圓形， 或者更多復雜的_結構元素_。

模糊/平滑


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
>>> lena = misc.lena()
>>> 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)
  
示例源碼 

消噪


向lena增加噪聲：
>>> from scipy import misc
>>> l = misc.lena()
>>> l = l[230:310, 210:350]
>>> 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)
  
示例源碼 
中值濾鏡：對直邊界效果更好(低曲率)：
>>> 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
其它局部非線性濾波器：維納濾波器(scipy.signal.wiener)等
非局部濾波器 
_總變差(TV)_消噪。找到新的圖像讓圖像的總變差(正態L1梯度的積分)變得最小，當接近測量圖像時：
>>> # from skimage.filter import tv_denoise
>>> from tv_denoise import tv_denoise
>>> tv_denoised = tv_denoise(noisy, weight=10)
>>> # More denoising (to the expense of fidelity to data)
>>> 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 
結構元素：
>>> 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]])
  
膨脹：最大化濾鏡:
>>> 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等等。

特征提取



邊緣檢測


合成數據：
>>> 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)
  
示例源碼 
canny濾鏡 
Canny濾鏡可以從skimage中獲取(文檔)，但是為了方便我們在這個教程中作為一個_單獨模塊_導入：
>>> #from skimage.filter import canny
>>> #or use module shipped with tutorial
>>> im += 0.1*np.random.random(im.shape)
>>> edges = canny(im, 1, 0.4, 0.2) # not enough smoothing
>>> edges = canny(im, 3, 0.3, 0.2) # better parameters
  
示例源碼 
需要調整幾個參數……過度擬合的風險

分割




基于_直方圖_的分割(沒有空間信息) 
  >>> 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
  
使用數學形態學來清理結果：
>>> # Remove small white regions
>>> open_img = ndimage.binary_opening(binary_img)
>>> # Remove small black hole
>>> close_img = ndimage.binary_closing(open_img)
  
示例源碼 
練習 
參看重建(reconstruction)操作(腐蝕+傳播(propagation))產生比開/閉操作更好的結果：
>>> 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.014678955078125
>>> np.abs(mask - reconstruct_final).mean()
0.0042572021484375
 練習 
檢查首次消噪步驟(中值濾波，總變差)如何更改直方圖，并且查看是否基于直方圖的分割更加精準了。


_基于圖像_的分割：使用空間信息 
  >>> 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, k=4, mode='arpack')
  >>> label_im = -np.ones(mask.shape)
  >>> label_im[mask] = labels

 

測量對象屬性：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?
  23
  >>> plt.imshow(label_im)        
  <matplotlib.image.AxesImage object at ...>

 
示例源碼 
計算每個區域的尺寸，均值等等：
>>> 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 ...>
 現在使用np.searchsorted重新分配標簽：
>>> labels = np.unique(label_im)
>>> label_im = np.searchsorted(labels, label_im)
  
示例源碼 
找到關注的封閉對象區域：³ 
>>> 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 ...>
  
示例源碼 
其它空間測量：ndiamge.center_of_mass,ndimage.maximum_position等等。
可以在分割應用限制范圍之外使用。
示例：塊平均(block mean)：
m scipy import misc
>>> l = misc.lena()
>>> sx, sy = l.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(l, labels=regions, index=np.arange(1,
...     regions.max() +1))
>>> block_mean.shape = (sx/4, sy/6)
 
您可能感興趣的文章