第 1 章 基本的圖像操作和處理
本章講解操作和處理圖像的基礎知識，將通過大量示例介紹處理圖像所需的 Python 工具包，並介紹用於讀取圖像、圖像轉換和縮放、計算導數、畫圖和保存結果等的基本工具。這些工具的使用將貫穿本書的剩余章節。
1.1 PIL：Python圖像處理類庫
PIL（Python Imaging Library Python，圖像處理類庫）提供了通用的圖像處理功能，以及大量有用的基本圖像操作，比如圖像縮放、裁剪、旋轉、顏色轉換等。PIL 是免費的，可以從 http://www.pythonware.com/products/pil/ 下載。
利用 PIL 中的函數，我們可以從大多數圖像格式的文件中讀取數據，然后寫入最常見的圖像格式文件中。PIL 中最重要的模塊為 Image。要讀取一幅圖像，可以使用：
from PIL import Image pil_im = Image.open('empire.jpg')上述代碼的返回值 pil_im 是一個 PIL 圖像對象。
圖像的顏色轉換可以使用 convert() 方法來實現。要讀取一幅圖像，並將其轉換成灰度圖像，只需要加上 convert('L')，如下所示：
pil_im = Image.open('empire.jpg').convert('L')在 PIL 文檔中有一些例子，參見 http://www.pythonware.com/library/pil/handbook/index.htm。這些例子的輸出結果如圖 1-1 所示。
圖 1-1：用 PIL 處理圖像的例子
1.1.1 轉換圖像格式
通過 save() 方法，PIL 可以將圖像保存成多種格式的文件。下面的例子從文件名列表（filelist）中讀取所有的圖像文件，並轉換成 JPEG 格式：
from PIL import Image import os for infile in filelist: outfile = os.path.splitext(infile)[0] + ".jpg" if infile != outfile: try: Image.open(infile).save(outfile) except IOError: print "cannot convert", infile PIL 的 open() 函數用於創建 PIL 圖像對象，save() 方法用於保存圖像到具有指定文件名的文件。除了后綴變為“.jpg”，上述代碼的新文件名和原文件名相同。PIL 是個足夠智能的類庫，可以根據文件擴展名來判定圖像的格式。PIL 函數會進行簡單的檢查，如果文件不是 JPEG 格式，會自動將其轉換成 JPEG 格式；如果轉換失敗，它會在控制台輸出一條報告失敗的消息。
本書會處理大量圖像列表。下面將創建一個包含文件夾中所有圖像文件的文件名列表。首先新建一個文件，命名為 imtools.py，來存儲一些經常使用的圖像操作，然后將下面的函數添加進去：
import os def get_imlist(path): """ 返回目錄中所有JPG 圖像的文件名列表""" return [os.path.join(path,f) for f in os.listdir(path) if f.endswith('.jpg')]現在，回到 PIL。
1.1.2 創建縮略圖
使用 PIL 可以很方便地創建圖像的縮略圖。thumbnail() 方法接受一個元組參數（該參數指定生成縮略圖的大小），然后將圖像轉換成符合元組參數指定大小的縮略圖。例如，創建最長邊為 128 像素的縮略圖，可以使用下列命令：
pil_im.thumbnail((128,128))1.1.3 復制和粘貼圖像區域
使用 crop() 方法可以從一幅圖像中裁剪指定區域：
box = (100,100,400,400) region = pil_im.crop(box)該區域使用四元組來指定。四元組的坐標依次是（左，上，右，下）。PIL 中指定坐標系的左上角坐標為（0，0）。我們可以旋轉上面代碼中獲取的區域，然后使用 paste() 方法將該區域放回去，具體實現如下：
region = region.transpose(Image.ROTATE_180) pil_im.paste(region,box)1.1.4 調整尺寸和旋轉
要調整一幅圖像的尺寸，我們可以調用 resize() 方法。該方法的參數是一個元組，用來指定新圖像的大小：
out = pil_im.resize((128,128))要旋轉一幅圖像，可以使用逆時針方式表示旋轉角度，然后調用 rotate() 方法：
out = pil_im.rotate(45)上述例子的輸出結果如圖 1-1 所示。最左端是原始圖像，然后是灰度圖像、粘貼有旋轉后裁剪圖像的原始圖像，最后是縮略圖。
1.2 Matplotlib
我們處理數學運算、繪制圖表，或者在圖像上繪制點、直線和曲線時，Matplotlib 是個很好的類庫，具有比 PIL 更強大的繪圖功能。Matplotlib 可以繪制出高質量的圖表，就像本書中的許多插圖一樣。Matplotlib 中的 PyLab 接口包含很多方便用戶創建圖像的函數。Matplotlib 是開源工具，可以從 http://matplotlib.sourceforge.net/ 免費下載。該鏈接中包含非常詳盡的使用說明和教程。下面的例子展示了本書中需要使用的大部分函數。
1.2.1 繪制圖像、點和線
盡管 Matplotlib 可以繪制出較好的條形圖、餅狀圖、散點圖等，但是對於大多數計算機視覺應用來說，僅僅需要用到幾個繪圖命令。最重要的是，我們想用點和線來表示一些事物，比如興趣點、對應點以及檢測出的物體。下面是用幾個點和一條線繪制圖像的例子：
from PIL import Image from pylab import * # 讀取圖像到數組中 im = array(Image.open('empire.jpg')) # 繪制圖像 imshow(im) # 一些點 x = [100,100,400,400] y = [200,500,200,500] # 使用紅色星狀標記繪制點 plot(x,y,'r*') # 繪制連接前兩個點的線 plot(x[:2],y[:2]) # 添加標題，顯示繪制的圖像 title('Plotting: "empire.jpg"') show()上面的代碼首先繪制出原始圖像，然后在 x 和 y 列表中給定點的 x 坐標和 y 坐標上繪制出紅色星狀標記點，最后在兩個列表表示的前兩個點之間繪制一條線段（默認為藍色）。該例子的繪制結果如圖 1-2 所示。show() 命令首先打開圖形用戶界面（GUI），然后新建一個圖像窗口。該圖形用戶界面會循環阻斷腳本，然后暫停，直到最后一個圖像窗口關閉。在每個腳本里，你只能調用一次 show() 命令，而且通常是在腳本的結尾調用。注意，在 PyLab 庫中，我們約定圖像的左上角為坐標原點。
圖像的坐標軸是一個很有用的調試工具；但是，如果你想繪制出較美觀的圖像，加上下列命令可以使坐標軸不顯示：
axis('off')上面的命令將繪制出如圖 1-2 右邊所示的圖像。
圖 1-2：Matplotlib 繪圖示例。帶有坐標軸和不帶坐標軸的包含點和一條線段的圖像
在繪圖時，有很多選項可以控制圖像的顏色和樣式。最有用的一些短命令如表 1-1、表 1-2 和表 1-3 所示。使用方法見下面的例子：
plot(x,y) # 默認為藍色實線 plot(x,y,'r*') # 紅色星狀標記 plot(x,y,'go-') # 帶有圓圈標記的綠線 plot(x,y,'ks:') # 帶有正方形標記的黑色虛線表1-1：用PyLab庫繪圖的基本顏色格式命令
| 顏色 |
|
|---|---|
|
|
藍色 |
|
|
綠色 |
|
|
紅色 |
|
|
青色 |
|
|
品紅 |
|
|
黃色 |
|
|
黑色 |
|
|
白色 |
表1-2：用PyLab庫繪圖的基本線型格式命令
| 線型 |
|
|---|---|
|
|
實線 |
|
|
虛線 |
|
|
點線 |
表1-3：用PyLab庫繪圖的基本繪制標記格式命令
| 標記 |
|
|---|---|
|
|
點 |
|
|
圓圈 |
|
|
正方形 |
|
|
星形 |
|
|
加號 |
|
|
叉號 |
1.2.2 圖像輪廓和直方圖
下面來看兩個特別的繪圖示例：圖像的輪廓和直方圖。繪制圖像的輪廓（或者其他二維函數的等輪廓線）在工作中非常有用。因為繪制輪廓需要對每個坐標 [x, y] 的像素值施加同一個閾值，所以首先需要將圖像灰度化：
from PIL import Image from pylab import * # 讀取圖像到數組中 im = array(Image.open('empire.jpg').convert('L')) # 新建一個圖像 figure() # 不使用顏色信息 gray() # 在原點的左上角顯示輪廓圖像 contour(im, origin='image') axis('equal') axis('off')像之前的例子一樣，這里用 PIL 的 convert() 方法將圖像轉換成灰度圖像。
圖像的直方圖用來表征該圖像像素值的分布情況。用一定數目的小區間（bin）來指定表征像素值的范圍，每個小區間會得到落入該小區間表示范圍的像素數目。該（灰度）圖像的直方圖可以使用 hist() 函數繪制：
figure() hist(im.flatten(),128) show()hist() 函數的第二個參數指定小區間的數目。需要注意的是，因為 hist() 只接受一維數組作為輸入，所以我們在繪制圖像直方圖之前，必須先對圖像進行壓平處理。flatten() 方法將任意數組按照行優先准則轉換成一維數組。圖 1-3 為等輪廓線和直方圖圖像。
圖 1-3：用 Matplotlib 繪制圖像等輪廓線和直方圖
1.2.3 交互式標注
有時用戶需要和某些應用交互，例如在一幅圖像中標記一些點，或者標注一些訓練數據。PyLab 庫中的 ginput() 函數就可以實現交互式標注。下面是一個簡短的例子：
from PIL import Image from pylab import * im = array(Image.open('empire.jpg')) imshow(im) print 'Please click 3 points' x = ginput(3) print 'you clicked:',x show()上面的腳本首先繪制一幅圖像，然后等待用戶在繪圖窗口的圖像區域點擊三次。程序將這些點擊的坐標 [x, y] 自動保存在 x 列表里。
1.3 NumPy
NumPy（http://www.scipy.org/NumPy/）是非常有名的 Python 科學計算工具包，其中包含了大量有用的思想，比如數組對象（用來表示向量、矩陣、圖像等）以及線性代數函數。NumPy 中的數組對象幾乎貫穿用於本書的所有例子中 1 數組對象可以幫助你實現數組中重要的操作，比如矩陣乘積、轉置、解方程系統、向量乘積和歸一化，這為圖像變形、對變化進行建模、圖像分類、圖像聚類等提供了基礎。
1PyLab 實際上包含 NumPy 的一些內容，如數組類型。這也是我們能夠在 1.2 節使用數組類型的原因。
NumPy 可以從 http://www.scipy.org/Download 免費下載，在線說明文檔（http://docs.scipy.org/doc/numpy/）包含了你可能遇到的大多數問題的答案。關於 NumPy 的更多內容，請參考開源書籍 [24]。
1.3.1 圖像數組表示
在先前的例子中，當載入圖像時，我們通過調用 array() 方法將圖像轉換成 NumPy 的數組對象，但當時並沒有進行詳細介紹。NumPy 中的數組對象是多維的，可以用來表示向量、矩陣和圖像。一個數組對象很像一個列表（或者是列表的列表），但是數組中所有的元素必須具有相同的數據類型。除非創建數組對象時指定數據類型，否則數據類型會按照數據的類型自動確定。
對於圖像數據，下面的例子闡述了這一點：
im = array(Image.open('empire.jpg')) print im.shape, im.dtype im = array(Image.open('empire.jpg').convert('L'),'f') print im.shape, im.dtype 控制台輸出結果如下所示：
(800, 569, 3) uint8 (800, 569) float32 每行的第一個元組表示圖像數組的大小（行、列、顏色通道），緊接着的字符串表示數組元素的數據類型。因為圖像通常被編碼成無符號八位整數（uint8），所以在第一種情況下，載入圖像並將其轉換到數組中，數組的數據類型為“uint8”。在第二種情況下，對圖像進行灰度化處理，並且在創建數組時使用額外的參數“f”；該參數將數據類型轉換為浮點型。關於更多數據類型選項，可以參考圖書 [24]。注意，由於灰度圖像沒有顏色信息，所以在形狀元組中，它只有兩個數值。
數組中的元素可以使用下標訪問。位於坐標 i、j，以及顏色通道 k 的像素值可以像下面這樣訪問：
value = im[i,j,k]多個數組元素可以使用數組切片方式訪問。切片方式返回的是以指定間隔下標訪問該數組的元素值。下面是有關灰度圖像的一些例子：
im[i,:] = im[j,:] # 將第 j 行的數值賦值給第 i 行 im[:,i] = 100 # 將第 i 列的所有數值設為100 im[:100,:50].sum() # 計算前100 行、前 50 列所有數值的和 im[50:100,50:100] # 50~100 行，50~100 列（不包括第 100 行和第 100 列） im[i].mean() # 第 i 行所有數值的平均值 im[:,-1] # 最后一列 im[-2,:] (or im[-2]) # 倒數第二行注意，示例僅僅使用一個下標訪問數組。如果僅使用一個下標，則該下標為行下標。注意，在最后幾個例子中，負數切片表示從最后一個元素逆向計數。我們將會頻繁地使用切片技術訪問像素值，這也是一個很重要的思想。
我們有很多操作和方法來處理數組對象。本書將在使用到的地方逐一介紹。你可以查閱在線文檔或者開源圖書 [24] 獲取更多信息。
1.3.2 灰度變換
將圖像讀入 NumPy 數組對象后，我們可以對它們執行任意數學操作。一個簡單的例子就是圖像的灰度變換。考慮任意函數 f，它將 0...255 區間（或者 0...1 區間）映射到自身（意思是說，輸出區間的范圍和輸入區間的范圍相同）。下面是關於灰度變換的一些例子：
from PIL import Image from numpy import * im = array(Image.open('empire.jpg').convert('L')) im2 = 255 - im # 對圖像進行反相處理 im3 = (100.0/255) * im + 100 # 將圖像像素值變換到100...200 區間 im4 = 255.0 * (im/255.0)**2 # 對圖像像素值求平方后得到的圖像第一個例子將灰度圖像進行反相處理；第二個例子將圖像的像素值變換到 100...200 區間；第三個例子對圖像使用二次函數變換，使較暗的像素值變得更小。圖 1-4 為所使用的變換函數圖像。圖 1-5 是輸出的圖像結果。你可以使用下面的命令查看圖像中的最小和最大像素值：
print int(im.min()), int(im.max())
圖 1-4：灰度變換示例。三個例子中所使用函數的圖像，其中虛線表示恆等變換
圖 1-5：灰度變換。對圖像應用圖 1-4 中的函數：f(x)=255-x 對圖像進行反相處理（左）；f(x)=(100/255)x+100 對圖像進行變換（中）；f(x)=255(x/255)2 對圖像做二次變換（右）
如果試着對上面例子查看最小值和最大值，可以得到下面的輸出結果：
2 255 0 253 100 200 0 255array() 變換的相反操作可以使用 PIL 的 fromarray() 函數完成：
pil_im = Image.fromarray(im)如果你通過一些操作將“uint8”數據類型轉換為其他數據類型，比如之前例子中的 im3 或者 im4，那么在創建 PIL 圖像之前，需要將數據類型轉換回來：
pil_im = Image.fromarray(uint8(im))如果你並不十分確定輸入數據的類型，安全起見，應該先轉換回來。注意，NumPy 總是將數組數據類型轉換成能夠表示數據的“最低”數據類型。對浮點數做乘積或除法操作會使整數類型的數組變成浮點類型。
1.3.3 圖像縮放
NumPy 的數組對象是我們處理圖像和數據的主要工具。想要對圖像進行縮放處理沒有現成簡單的方法。我們可以使用之前 PIL 對圖像對象轉換的操作，寫一個簡單的用於圖像縮放的函數。把下面的函數添加到 imtool.py 文件里：
def imresize(im,sz): """ 使用PIL 對象重新定義圖像數組的大小""" pil_im = Image.fromarray(uint8(im)) return array(pil_im.resize(sz))我們將會在接下來的內容中使用這個函數。
1.3.4 直方圖均衡化
圖像灰度變換中一個非常有用的例子就是直方圖均衡化。直方圖均衡化是指將一幅圖像的灰度直方圖變平，使變換后的圖像中每個灰度值的分布概率都相同。在對圖像做進一步處理之前，直方圖均衡化通常是對圖像灰度值進行歸一化的一個非常好的方法，並且可以增強圖像的對比度。
在這種情況下，直方圖均衡化的變換函數是圖像中像素值的累積分布函數（cumulative distribution function，簡寫為 cdf，將像素值的范圍映射到目標范圍的歸一化操作）。
下面的函數是直方圖均衡化的具體實現。將這個函數添加到 imtool.py 里：
def histeq(im,nbr_bins=256): """ 對一幅灰度圖像進行直方圖均衡化""" # 計算圖像的直方圖 imhist,bins = histogram(im.flatten(),nbr_bins,normed=True) cdf = imhist.cumsum() # cumulative distribution function cdf = 255 * cdf / cdf[-1] # 歸一化 # 使用累積分布函數的線性插值，計算新的像素值 im2 = interp(im.flatten(),bins[:-1],cdf) return im2.reshape(im.shape), cdf 該函數有兩個輸入參數，一個是灰度圖像，一個是直方圖中使用小區間的數目。函數返回直方圖均衡化后的圖像，以及用來做像素值映射的累積分布函數。注意，函數中使用到累積分布函數的最后一個元素（下標為 -1），目的是將其歸一化到 0...1 范圍。你可以像下面這樣使用該函數：
from PIL import Image from numpy import * im = array(Image.open('AquaTermi_lowcontrast.jpg').convert('L')) im2,cdf = imtools.histeq(im)圖 1-6 和圖 1-7 為上面直方圖均衡化例子的結果。上面一行顯示的分別是直方圖均衡化之前和之后的灰度直方圖，以及累積概率分布函數映射圖像。可以看到，直方圖均衡化后圖像的對比度增強了，原先圖像灰色區域的細節變得清晰。
圖 1-6：直方圖均衡化示例。左側為原始圖像和直方圖，中間圖為灰度變換函數，右側為直方圖均衡化后的圖像和相應直方圖
圖 1-7：直方圖均衡化示例。左側為原始圖像和直方圖，中間圖為灰度變換函數，右側為直方圖均衡化后的圖像和相應直方圖
1.3.5 圖像平均
圖像平均操作是減少圖像噪聲的一種簡單方式，通常用於藝術特效。我們可以簡單地從圖像列表中計算出一幅平均圖像。假設所有的圖像具有相同的大小，我們可以將這些圖像簡單地相加，然后除以圖像的數目，來計算平均圖像。下面的函數可以用於計算平均圖像，將其添加到 imtool.py 文件里：
def compute_average(imlist): """ 計算圖像列表的平均圖像""" # 打開第一幅圖像，將其存儲在浮點型數組中 averageim = array(Image.open(imlist[0]), 'f') for imname in imlist[1:]: try: averageim += array(Image.open(imname)) except: print imname + '...skipped' averageim /= len(imlist) # 返回uint8 類型的平均圖像 return array(averageim, 'uint8')該函數包括一些基本的異常處理技巧，可以自動跳過不能打開的圖像。我們還可以使用 mean() 函數計算平均圖像。mean() 函數需要將所有的圖像堆積到一個數組中；也就是說，如果有很多圖像，該處理方式需要占用很多內存。我們將會在下一節中使用該函數。
1.3.6 圖像的主成分分析（PCA）
PCA（Principal Component Analysis，主成分分析）是一個非常有用的降維技巧。它可以在使用盡可能少維數的前提下，盡量多地保持訓練數據的信息，在此意義上是一個最佳技巧。即使是一幅 100×100 像素的小灰度圖像，也有 10 000 維，可以看成 10 000 維空間中的一個點。一兆像素的圖像具有百萬維。由於圖像具有很高的維數，在許多計算機視覺應用中，我們經常使用降維操作。PCA 產生的投影矩陣可以被視為將原始坐標變換到現有的坐標系，坐標系中的各個坐標按照重要性遞減排列。
為了對圖像數據進行 PCA 變換，圖像需要轉換成一維向量表示。我們可以使用 NumPy 類庫中的 flatten() 方法進行變換。
將變平的圖像堆積起來，我們可以得到一個矩陣，矩陣的一行表示一幅圖像。在計算主方向之前，所有的行圖像按照平均圖像進行了中心化。我們通常使用 SVD（Singular Value Decomposition，奇異值分解）方法來計算主成分；但當矩陣的維數很大時，SVD 的計算非常慢，所以此時通常不使用 SVD 分解。下面就是 PCA 操作的代碼：
from PIL import Image from numpy import * def pca(X): """ 主成分分析： 輸入：矩陣X ，其中該矩陣中存儲訓練數據，每一行為一條訓練數據 返回：投影矩陣（按照維度的重要性排序）、方差和均值""" # 獲取維數 num_data,dim = X.shape # 數據中心化 mean_X = X.mean(axis=0) X = X - mean_X if dim>num_data: # PCA- 使用緊致技巧 M = dot(X,X.T) # 協方差矩陣 e,EV = linalg.eigh(M) # 特征值和特征向量 tmp = dot(X.T,EV).T # 這就是緊致技巧 V = tmp[::-1] # 由於最后的特征向量是我們所需要的，所以需要將其逆轉 S = sqrt(e)[::-1] # 由於特征值是按照遞增順序排列的，所以需要將其逆轉 for i in range(V.shape[1]): V[:,i] /= S else: # PCA- 使用SVD 方法 U,S,V = linalg.svd(X) V = V[:num_data] # 僅僅返回前nun_data 維的數據才合理 # 返回投影矩陣、方差和均值 return V,S,mean_X 該函數首先通過減去每一維的均值將數據中心化，然后計算協方差矩陣對應最大特征值的特征向量，此時可以使用簡明的技巧或者 SVD 分解。這里我們使用了 range() 函數，該函數的輸入參數為一個整數 n，函數返回整數 0...(n-1) 的一個列表。你也可以使用 arange() 函數來返回一個數組，或者使用 xrange() 函數返回一個產生器（可能會提升速度）。我們在本書中貫穿使用 range() 函數。
如果數據個數小於向量的維數，我們不用 SVD 分解，而是計算維數更小的協方差矩陣 XXT 的特征向量。通過僅計算對應前 k（k 是降維后的維數）最大特征值的特征向量，可以使上面的 PCA 操作更快。由於篇幅所限，有興趣的讀者可以自行探索。矩陣 V 的每行向量都是正交的，並且包含了訓練數據方差依次減少的坐標方向。
我們接下來對字體圖像進行 PCA 變換。fontimages.zip 文件包含采用不同字體的字符 a 的縮略圖。所有的 2359 種字體可以免費下載 2。假定這些圖像的名稱保存在列表 imlist 中，跟之前的代碼一起保存傳在 pca.py 文件中，我們可以使用下面的腳本計算圖像的主成分：
2免費字體圖像庫由 Martin Solli 收集並上傳（http://webstaff.itn.liu.se/~marso/）。
from PIL import Image from numpy import * from pylab import * import pca im = array(Image.open(imlist[0])) # 打開一幅圖像，獲取其大小 m,n = im.shape[0:2] # 獲取圖像的大小 imnbr = len(imlist) # 獲取圖像的數目 # 創建矩陣，保存所有壓平后的圖像數據 immatrix = array([array(Image.open(im)).flatten() for im in imlist],'f') # 執行 PCA 操作 V,S,immean = pca.pca(immatrix) # 顯示一些圖像（均值圖像和前 7 個模式） figure() gray() subplot(2,4,1) imshow(immean.reshape(m,n)) for i in range(7): subplot(2,4,i+2) imshow(V[i].reshape(m,n)) show()注意，圖像需要從一維表示重新轉換成二維圖像；可以使用 reshape() 函數。如圖 1-8 所示，運行該例子會在一個繪圖窗口中顯示 8 個圖像。這里我們使用了 PyLab 庫的 subplot() 函數在一個窗口中放置多個圖像。
圖 1-8：平均圖像（左上）和前 7 個模式（具有最大方差的方向模式）
1.3.7 使用pickle模塊
如果想要保存一些結果或者數據以方便后續使用，Python 中的 pickle 模塊非常有用。pickle 模塊可以接受幾乎所有的 Python 對象，並且將其轉換成字符串表示，該過程叫做封裝（pickling）。從字符串表示中重構該對象，稱為拆封（unpickling）。這些字符串表示可以方便地存儲和傳輸。
我們來看一個例子。假設想要保存上一節字體圖像的平均圖像和主成分，可以這樣來完成：
# 保存均值和主成分數據 f = open('font_pca_modes.pkl', 'wb') pickle.dump(immean,f) pickle.dump(V,f) f.close()在上述例子中，許多對象可以保存到同一個文件中。pickle 模塊中有很多不同的協議可以生成 .pkl 文件；如果不確定的話，最好以二進制文件的形式讀取和寫入。在其他 Python 會話中載入數據，只需要如下使用 load() 方法：
# 載入均值和主成分數據 f = open('font_pca_modes.pkl', 'rb') immean = pickle.load(f) V = pickle.load(f) f.close()注意，載入對象的順序必須和先前保存的一樣。Python 中有個用 C 語言寫的優化版本，叫做 cpickle 模塊，該模塊和標准 pickle 模塊完全兼容。關於 pickle 模塊的更多內容，參見 pickle 模塊文檔頁 http://docs.python.org/library/pickle.html。
在本書接下來的章節中，我們將使用 with 語句處理文件的讀寫操作。這是 Python 2.5 引入的思想，可以自動打開和關閉文件（即使在文件打開時發生錯誤）。下面的例子使用 with() 來實現保存和載入操作：
# 打開文件並保存 with open('font_pca_modes.pkl', 'wb') as f: pickle.dump(immean,f) pickle.dump(V,f)和
# 打開文件並載入 with open('font_pca_modes.pkl', 'rb') as f: immean = pickle.load(f) V = pickle.load(f)上面的例子乍看起來可能很奇怪，但 with() 確實是個很有用的思想。如果你不喜歡它，可以使用之前的 open 和 close函數。
作為 pickle 的一種替代方式，NumPy 具有讀寫文本文件的簡單函數。如果數據中不包含復雜的數據結構，比如在一幅圖像上點擊的點列表，NumPy 的讀寫函數會很有用。保存一個數組 x 到文件中，可以使用：
savetxt('test.txt',x,'%i')最后一個參數表示應該使用整數格式。類似地，讀取可以使用：
x = loadtxt('test.txt')你可以從在線文檔 http://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html 了解更多內容。
最后，NumPy 有專門用於保存和載入數組的函數。你可以在上面的在線文檔里查看關於 save() 和 load() 的更多內容。
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 Image from numpy import * from scipy.ndimage import filters im = 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 Image from numpy import * from scipy.ndimage import filters im = 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<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。
-
讀寫.mat文件
如果你有一些數據，或者在網上下載到一些有趣的數據集，這些數據以 Matlab 的 .mat 文件格式存儲，那么可以使用
scipy.io模塊進行讀取。data = scipy.io.loadmat('test.mat')上面代碼中，
data對象包含一個字典，字典中的鍵對應於保存在原始 .mat 文件中的變量名。由於這些變量是數組格式的，因此可以很方便地保存到 .mat 文件中。你僅需創建一個字典（其中要包含你想要保存的所有變量），然后使用savemat()函數：data = {} data['x'] = x scipy.io.savemat('test.mat',data)因為上面的腳本保存的是數組 x，所以當讀入到 Matlab 中時，變量的名字仍為 x。關於
scipy.io模塊的更多內容，請參見在線文檔 http://docs.scipy.org/doc/scipy/reference/io.html。 -
以圖像形式保存數組
因為我們需要對圖像進行操作，並且需要使用數組對象來做運算，所以將數組直接保存為圖像文件 4 非常有用。本書中的很多圖像都是這樣的創建的。
imsave()函數可以從scipy.misc模塊中載入。要將數組im保存到文件中，可以使用下面的命令：from scipy.misc import imsave imsave('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_init Px = im # 對偶域的x 分量 Py = im # 對偶域的y 分量 error = 1 while (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*