Python影像處理庫(2)

標籤：

1.4　SciPy

SciPy（http://scipy.org/） 是建立在 NumPy 基礎上，用於數值運算的開源工具包。SciPy 提供很多高效的操作，可以實現數值積分、最佳化、統計、訊號處理，以及對我們來說最重要的影像處理功能。接下來，本節會介紹 SciPy 中大量有用的模組。SciPy 是個開源工具包，可以從http://scipy.org/Download 下載。

1.4.1　映像模糊

映像的高斯模糊是非常經典的映像卷積例子。本質上，映像模糊就是將（灰階）映像 I 和一個高斯核進行卷積操作：

Iσ = I*Gσ

其中 * 表示卷積操作；Gσ 是標準差為 σ 的二維高斯核，定義為 :

高斯模糊通常是其他影像處理操作的一部分，比像插值操作、興趣點計算以及很多其他應用。

SciPy 有用來做濾波操作的 scipy.ndimage.filters 模組。該模組使用快速一維分離的方式來計算卷積。你可以像下面這樣來使用它：

from PIL import Imagefrom numpy import *from scipy.ndimage import filtersim = array(Image.open(‘empire.jpg‘).convert(‘L‘))im2 = filters.gaussian_filter(im,5)

上面 guassian_filter() 函數的最後一個參數表示標準差。

圖 1-9 顯示了隨著 σ 的增加，一幅映像被模糊的程度。σ 越大，處理後的映像細節丟失越多。如果打算模糊一幅彩色映像，只需簡單地對每一個色彩通道進行高斯模糊：

im = array(Image.open(‘empire.jpg‘))im2 = zeros(im.shape)for i in range(3):  im2[:,:,i] = filters.gaussian_filter(im[:,:,i],5)im2 = uint8(im2)

在上面的指令碼中，最後並不總是需要將映像轉換成 uint8 格式，這裡只是將像素值用八位來表示。我們也可以使用：

im2 = array(im2,‘uint8‘)

來完成轉換。

關於該模組更多的內容以及不同參數的選擇，請查看http://docs.scipy.org/doc/scipy/reference/ndimage.html 上 SciPy 文檔中的 scipy.ndimage部分。

圖 1-9：使用 scipy.ndimage.filters 模組進行高斯模糊：（a）原始灰階映像；（b）使用 σ=2 的高斯濾波器；（c）使用 σ=5 的高斯濾波器；（d）使用 σ=10 的高斯濾波器

1.4.2　映像導數

整本書中可以看到，在很多應用中映像強度的變化情況是非常重要的資訊。強度的變化可以用灰階映像 I（對於彩色映像，通常對每個色彩通道分別計算導數）的 x 和 y 方嚮導數 Ix 和 Iy 進行描述。

映像的梯度向量為∇I = [Ix, Iy]T。梯度有兩個重要的屬性，一是梯度的大小：

它描述了映像強度變化的強弱，一是梯度的角度：

α=arctan2(Iy, Ix)

描述了映像中在每個點（像素）上強度變化最大的方向。NumPy 中的 arctan2() 函數返回弧度表示的有符號角度，角度的變化區間為 -π...π。

我們可以用離散近似的方式來計算映像的導數。映像導數大多數可以通過卷積簡單地實現：

Ix=I*Dx 和 Iy=I*Dy

對於 Dx 和 Dy，通常選擇 Prewitt 濾波器：

 和 

或者 Sobel 濾波器：

 和 

這些導數濾波器可以使用 scipy.ndimage.filters 模組的標準卷積操作來簡單地實現，例如：

from PIL import Imagefrom numpy import *from scipy.ndimage import filtersim = array(Image.open(‘empire.jpg‘).convert(‘L‘))# Sobel 導數濾波器imx = zeros(im.shape)filters.sobel(im,1,imx)imy = zeros(im.shape)filters.sobel(im,0,imy)magnitude = sqrt(imx**2+imy**2)

上面的指令碼使用 Sobel 濾波器來計算 x 和 y 的方嚮導數，以及梯度大小。sobel() 函數的第二個參數表示選擇 x 或者 y 方嚮導數，第三個參數儲存輸出的變數。圖 1-10 顯示了用 Sobel 濾波器計算出的導數映像。在兩個導數映像中，正導數顯示為亮的像素，負導數顯示為暗的像素。灰色地區表示導數的值接近於零。

圖 1-10：使用 Sobel 導數濾波器計算導數映像：（a）原始灰階映像；（b）x 導數映像；（c）y導數映像；（d）梯度大小映像

上述計算映像導數的方法有一些缺陷：在該方法中，濾波器的尺度需要隨著映像解析度的變化而變化。為了在映像雜訊方面更穩健，以及在任意尺度上計算導數，我們可以使用高斯導數濾波器：

Ix=I*Gσx 和 Iy=I*Gσy

其中，Gσx和 Gσy 表示 Gσ 在 x 和 y 方向上的導數，Gσ 為標準差為 σ 的高斯函數。

我們之前用於模糊的 filters.gaussian_filter() 函數可以接受額外的參數，用來計算高斯導數。可以簡單地按照下面的方式來處理：

sigma = 5 # 標準差imx = zeros(im.shape)filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)imy = zeros(im.shape)filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)

該函數的第三個參數指定對每個方向計算哪種類型的導數，第二個參數為使用的標準差。你可以查看相應文檔瞭解詳情。圖 1-11 顯示了不同尺度下的導數映像和梯度大小。你可以和圖 1-9 中做相同尺度模糊的映像做比較。

圖 1-11：使用高斯導數計算映像導數：x 導數映像（上），y 導數映像（中），以及梯度大小映像（下）；（a）為原始灰階映像，（b）為使用 σ=2 的高斯導數濾波器處理後的映像，（c）為使 用 σ=5 的高斯導數濾波器處理後的映像，（d）為使用 σ=10 的高斯導數濾波器處理後的映像

1.4.3　形態學：對象計數

形態學（或數學形態學）是度量和分析基本形狀的影像處理方法的基本架構與集合。形態學通常用於處理二值映像，但是也能夠用於灰階映像。二值映像是指映像的每個像素只能取兩個值，通常是 0 和 1。二值映像通常是，在計算物體的數目，或者度量其大小時，對一幅映像進行閾值化後的結果。你可以從 http://en.wikipedia.org/wiki/Mathematical_morphology 大體瞭解形態學及其處理映像的方式。

scipy.ndimage 中的 morphology 模組可以實現形態學操作。你可以使用 scipy.ndimage 中的measurements 模組來實現二值映像的計數和度量功能。下面通過一個簡單的例子介紹如何使用它們。

考慮在圖 1-12a3 裡的二值映像，計算該映像中的對象個數可以通過下面的指令碼實現：

3這個映像實際上是映像“分割”後的結果。如果你想知道該映像是如何建立的，可以查看 9.3 節。

from scipy.ndimage import measurements,morphology# 載入映像，然後使用閾值化操作，以保證處理的映像為二值映像im = array(Image.open(‘houses.png‘).convert(‘L‘))im = 1*(im<</span>128)labels, nbr_objects = measurements.label(im)print "Number of objects:", nbr_objects

上面的指令碼首先載入該映像，通過閾值化方式來確保該映像是二值映像。通過和 1 相乘，指令碼將布爾數群組轉換成二進位表示。然後，我們使用 label() 函數尋找單個的物體，並且按照它們屬於哪個對象將整數標籤給像素賦值。圖 1-12b 是 labels 數組的映像。映像的灰階值表示對象的標籤。可以看到，在一些對象之間有一些小的串連。進行二進位開（binary open）操作，我們可以將其移除：

# 形態學開操作更好地分離各個對象im_open = morphology.binary_opening(im,ones((9,5)),iterations=2)labels_open, nbr_objects_open = measurements.label(im_open)print "Number of objects:", nbr_objects_open

binary_opening() 函數的第二個參數指定一個數組結構元素。該數組表示以一個像素為中心時，使用哪些相鄰像素。在這種情況下，我們在 y 方向上使用 9 個像素（上面 4 個像素、像素本身、下面 4 個像素），在 x 方向上使用 5 個像素。你可以指定任意數組為結構元素，數組中的非零元素決定使用哪些相鄰像素。參數 iterations 決定執行該操作的次數。你可以嘗試使用不同的迭代次數 iterations 值，看一下對象的數目如何變化。你可以在圖 1-12c 與圖 1-12d 中查看經過開操作後的映像，以及相應的標籤映像。正如你想象的一樣，binary_closing() 函數實現相反的操作。我們將該函數和在 morphology 和 measurements 模組中的其他函數的用法留作練習。你可以從 scipy.ndimage 模組文檔 http://docs.scipy.org/doc/scipy/reference/ndimage.html 中瞭解關於這些函數的更多知識。

圖 1-12：形態學樣本。使用二值開操作將對象分開，然後計算物體的數目：（a）為原始二值映像；（b）為對應原始映像的標籤映像，其中灰階值表示物體的標籤；（c）為使用開操作後的二值映像；（d）為開操作後映像的標籤映像

1.4.4　一些有用的SciPy模組

SciPy 中包含一些用於輸入和輸出的實用模組。下面介紹其中兩個模組：io 和 misc。

1. 讀寫.mat檔案

   如果你有一些資料，或者在網上下載到一些有趣的資料集，這些資料以 Matlab 的 .mat 檔案格式儲存，那麼可以使用 scipy.io 模組進行讀取。

   data = scipy.io.loadmat(‘test.mat‘)

   上面代碼中，data 對象包含一個字典，字典中的鍵對應於儲存在原始 .mat 檔案中的變數名。由於這些變數是數組格式的，因此可以很方便地儲存到 .mat 檔案中。你僅需建立一個字典（其中要包含你想要儲存的所有變數），然後使用 savemat() 函數：

   data = {}data[‘x‘] = xscipy.io.savemat(‘test.mat‘,data)

   因為上面的指令碼儲存的是數組 x，所以當讀入到 Matlab 中時，變數的名字仍為 x。關於scipy.io 模組的更多內容，請參見線上文檔http://docs.scipy.org/doc/scipy/reference/io.html。

2. 以映像形式儲存數組

   因為我們需要對映像進行操作，並且需要使用數組對象來做運算，所以將數組直接儲存為影像檔 4 非常有用。本書中的很多映像都是這樣的建立的。

   imsave() 函數可以從 scipy.misc 模組中載入。要將數組 im 儲存到檔案中，可以使用下面的命令：

   from scipy.misc import imsaveimsave(‘test.jpg‘,im)

   scipy.misc 模組同樣包含了著名的 Lena 測試映像：

   lena = scipy.misc.lena()

   該指令碼返回一個 512×512 的灰階映像數組。

4所有 Pylab 圖均可儲存為多種映像格式，方法是點擊映像視窗中的“儲存”按鈕。

1.5　進階樣本：映像去噪

我們通過一個非常實用的例子——映像的去噪——來結束本章。映像去噪是在去除映像雜訊的同時，儘可能地保留映像細節和結構的處理技術。我們這裡使用 ROF（Rudin-Osher-Fatemi）去噪模型。該模型最早出現在文獻 [28] 中。映像去噪對於很多應用來說都非常重要；這些應用範圍很廣，小到讓你的假期照片看起來更漂亮，大到提高衛星映像的品質。ROF 模型具有很好的性質：使處理後的映像更平滑，同時保持映像邊緣和結構資訊。

ROF 模型的數學基礎和處理技巧非常高深，不在本書講述範圍之內。在講述如何基於 Chambolle 提出的演算法 [5] 實現 ROF 求解器之前，本書首先簡要介紹一下 ROF 模型。

一幅（灰階）映像 I 的全變差（Total Variation，TV）定義為梯度範數之和。在連續表示的情況下，全變差表示為：

　　　　　　　　　　　　（1.1）

在離散表示的情況下，全變差表示為：

其中，上面的式子是在所有映像座標 x=[x, y] 上取和。

在 Chambolle 提出的 ROF 模型裡，目標函數為尋找降噪後的映像 U，使下式最小：

其中範數 ||I-U|| 是去噪後映像 U 和原始映像 I 差異的度量。也就是說，本質上該模型使去噪後的映像像素值“平坦”變化，但是在映像地區的邊緣上，允許去噪後的映像像素值“跳躍”變化。

按照論文 [5] 中的演算法，我們可以按照下面的代碼實現 ROF 模型去噪：

from numpy import *def denoise(im,U_init,tolerance=0.1,tau=0.125,tv_weight=100):  """ 使用A. Chambolle（2005）在公式（11）中的計算步驟實現Rudin-Osher-Fatemi（ROF）去噪模型    輸入：含有雜訊的輸入映像（灰階映像）、U 的初始值、TV 正則項權值、步長、停業條件    輸出：去噪和去除紋理後的映像、紋理殘留"""m,n = im.shape # 雜訊映像的大小# 初始化U = U_initPx = im # 對偶域的x 分量Py = im # 對偶域的y 分量error = 1while (error > tolerance):  Uold = U  # 原始變數的梯度  GradUx = roll(U,-1,axis=1)-U # 變數U 梯度的x 分量  GradUy = roll(U,-1,axis=0)-U # 變數U 梯度的y 分量  # 更新對偶變數  PxNew = Px + (tau/tv_weight)*GradUx  PyNew = Py + (tau/tv_weight)*GradUy  NormNew = maximum(1,sqrt(PxNew**2+PyNew**2))  Px = PxNew/NormNew # 更新x 分量（對偶）  Py = PyNew/NormNew # 更新y 分量（對偶）  # 更新原始變數  RxPx = roll(Px,1,axis=1) # 對x 分量進行向右x 軸平移  RyPy = roll(Py,1,axis=0) # 對y 分量進行向右y 軸平移  DivP = (Px-RxPx)+(Py-RyPy) # 對偶域的散度  U = im + tv_weight*DivP # 更新原始變數  # 更新誤差  error = linalg.norm(U-Uold)/sqrt(n*m);return U,im-U # 去噪後的映像和紋理殘餘

在這個例子中，我們使用了 roll() 函數。顧名思義，在一個座標軸上，它迴圈“滾動”數組中的元素值。該函數可以非常方便地計算鄰域元素的差異，比如這裡的導數。我們還使用了linalg.norm() 函數，該函數可以衡量兩個數組間（這個例子中是指映像矩陣 U和 Uold）的差異。我們將這個 denoise() 函數儲存到 rof.py 檔案中。

下面使用一個合成的雜訊映像樣本來說明如何使用該函數：

from numpy import *from numpy import randomfrom scipy.ndimage import filtersimport rof# 使用雜訊建立合成映像im = zeros((500,500))im[100:400,100:400] = 128im[200:300,200:300] = 255im = im + 30*random.standard_normal((500,500))U,T = rof.denoise(im,im)G = filters.gaussian_filter(im,10)# 儲存產生結果from scipy.misc import imsaveimsave(‘synth_rof.pdf‘,U)imsave(‘synth_gaussian.pdf‘,G)

原始映像和映像的去噪結果 1-13 所示。正如你所看到的，ROF 演算法去噪後的映像很好地保留了映像的邊緣資訊。

圖 1-13：使用 ROF 模型對合成映像去噪：（a）為原始雜訊映像；（b）為經過高斯模糊的映像（σ=10）；（c）為經過 ROF 模型去噪後的映像

下面看一下在實際映像中使用 ROF 模型去噪的效果：

from PIL import Imagefrom pylab import *import rofim = array(Image.open(‘empire.jpg‘).convert(‘L‘))U,T = rof.denoise(im,im)figure()gray()imshow(U)axis(‘equal‘)axis(‘off‘)show()

經過 ROF 去噪後的映像 1-14c 所示。為了方便比較，該圖中同樣顯示了模糊後的映像。可以看到，ROF 去噪後的映像保留了邊緣和映像的結構資訊，同時模糊了“雜訊”。

圖 1-14：使用 ROF 模型對灰階映像去噪：（a）為原始雜訊映像；（b）為經過高斯模糊的映像（σ=5）；（c）為經過 ROF 模型去噪後的映像

Python影像處理庫(2)