scipy.ndimage是一個處理多元圖像的函數庫，其中又包括以下幾個子產品：
- ·filters：圖像濾波器。
- ·fourier：傅立葉變換。
- ·interpolation：圖像的插值、旋轉以及仿射變換等。
- ·measurements：圖像相關資訊的測量。
- ·morphology：形态學圖像處理。
更強大的圖像處理庫
scipy.ndimage隻提供了一些基礎的圖像處理功能，下面是一些更強大的圖像處理庫：
·OpenCV：它是使用C/C++開發的計算機視覺庫，本書将用一整章的篇幅介紹OpenCV提供的Python調用接口的用法。
·SimpleCV：對多個計算機視覺庫進行包裝，提供了一套更友善、統一的Python調用接口。
·scikit-image：使用Python開發的圖像處理庫，高速運算部分多采用Cython編寫。
·Mahotas：采用Python和C++開發的圖像處理庫。
一、形态學圖像處理
本節介紹如何使用morphology子產品實作二值圖像處理。二值圖像中每個像素的顔色隻有兩種：黑色和白色。在NumPy中可以用二維布爾數組表示：False表示黑色，True表示白色。也可以用無符号單位元組整型(uint8)數組表示：0表示黑色，非0表示白色。
下面的兩個函數用于顯示形态學圖像處理的結果：
import numpy as np
def expand_image(img, value, out=None, size = 10):
if out is None:
w, h = img.shape
out = np.zeros((w*size, h*size),dtype=np.uint8)
tmp = np.repeat(np.repeat(img,size,0),size,1)
out[:,:] = np.where(tmp, value, out)
out[::size,:] = 0
out[:,::size] = 0
return out
def show_image(*imgs):
for idx, img in enumerate(imgs, 1):
ax = pl.subplot(1, len(imgs), idx)
pl.imshow(img, cmap="gray")
ax.set_axis_off( )
pl.subplots_adjust(0.02, 0, 0.98, 1, 0.02, 0) 1．膨脹和腐蝕
二值圖像最基本的形态學運算是膨脹和腐蝕。膨脹運算是将與某物體（白色區域）接觸的所有背景像素（黑色區域）合并到該物體中的過程。簡單來說，就是對于原始圖像中的每個白色像素進行處理，将其周圍的黑色像素都設定為白色像素。這裡的“周圍”是一個模糊概念，在實際運算時，需要明确給出“周圍”的定義。圖1是膨脹運算的一個例子，其中左圖是原始圖像，中間的圖是四連通定義的“周圍”的膨脹效果，右圖是八連通定義的“周圍”的膨脹效果。圖中用灰色方塊表示由膨脹處理添加進物體的像素。
from scipy.ndimage import morphology
def dilation_demo(a, structure=None):
b = morphology.binary_dilation(a, structure)
img = expand_image(a, 255)
return expand_image(np.logical_xor(a,b), 150, out=img)
a = pl.imread("scipy_morphology_demo.png")[:,:,0].astype(np.uint8)
img1 = expand_image(a, 255)
img2 = dilation_demo(a)
img3 = dilation_demo(a, [[1,1,1],[1,1,1],[1,1,1]])
show_image(img1, img2, img3) 
圖1 四連通和八連通的膨脹運算
四連通包括上下左右4個像素，而八連通則還包括4個斜線方向上的鄰接像素。它們都可以使用下面的正方形矩陣定義，其中正中心的元素表示目前要進行運算的像素，而其周圍的1和0表示對應位置的像素是否算作其“周圍”像素。這種矩陣描述了周圍像素和目前像素之間的關系，被稱作結構元素(structuring element)。
四連通八連通
0 1 0 1 1 1
1 1 1 1 1 1
0 1 0 1 1 1
假設數組a是一個表示二值圖像的數組，可以用如下語句對其進行膨脹運算：
binary_dilation(a)
binary_dilation( )預設使用四連通進行膨脹運算，通過structure參數可以指定其他的結構元素。下面是進行八連通膨脹運算的語句：
binary_dilation(a, structure=[[1,1,1],[1,1,1],[1,1,1]])
通過設定不同的結構元素，能夠實作各種不同的效果，圖2顯示了三種不同結構元素的膨脹效果。圖中的結構元素分别為：
圖2 不同結構元素的膨脹效果
左中右
0 0 0 0 1 0 0 1 0
1 1 1 0 1 0 0 1 0
0 0 0 0 1 0 0 0 0
img4 = dilation_demo(a, [[0,0,0],[1,1,1],[0,0,0]])
img5 = dilation_demo(a, [[0,1,0],[0,1,0],[0,1,0]])
img6 = dilation_demo(a, [[0,1,0],[0,1,0],[0,0,0]])
show_image(img4, img5, img6) binary_erosion()的腐蝕運算正好和膨脹相反，它将“周圍”有黑色像素的白色像素設定為黑色。圖3是四連通和八連通腐蝕的效果，圖中用灰色方塊表示被腐蝕的像素。
圖3 四連通和八連通的腐蝕運算
def erosion_demo(a, structure=None):
b = morphology.binary_erosion(a, structure)
img = expand_image(a, 255)
return expand_image(np.logical_xor(a,b), 100, out=img)
img1 = expand_image(a, 255)
img2 = erosion_demo(a)
img3 = erosion_demo(a, [[1,1,1],[1,1,1],[1,1,1]])
show_image(img1, img2, img3) 2．Hit和Miss
Hit和Miss是二值形态學圖像進行中最基本的運算，因為幾乎所有其他的運算都可以用Hit和Miss的組合推演出來。它對圖像中的每個像素周圍的像素進行模式判斷，如果周圍像素的黑白模式符合指定的模式，将此像素設為白色，否則設定為黑色。因為它需要同時對白色和黑色像素進行判斷，是以需要指定兩個結構元素。進行Hit和Miss運算的binary_hit_or_miss()的調用形式如下：
binary_hit_or_miss(input, structure1=None, structure2=None, ...)
其中structure1參數指定白色像素的結構元素，而structure2參數則指定黑色像素的結構元素。圖4是binary_hit_or_miss()的運算結果。其中左圖為原始圖像，中圖為使用下面兩個結構元素進行運算的結果：
圖4 Hit和Miss運算
白色結構元素 黑色結構元素
0 0 0 1 0 0
0 1 0 0 0 0
1 1 1 0 0 0
在這兩個結構元素中，0表示不關心其對應位置的像素的顔色，1表示其對應位置的像素必須為結構元素所表示的顔色。是以通過這兩個結構元素可以找到“下方三個像素為白色，并且左上像素為黑色的白色像素”。
與右圖對應的結構元素如下，通過它可以找到“下方三個像素為白色、左上像素為黑色的黑色像素”。
白色結構元素 黑色結構元素
0 0 0 1 0 0
0 0 0 0 1 0
1 1 1 0 0 0
def hitmiss_demo(a, structure1, structure2):
b = morphology.binary_hit_or_miss(a, structure1, structure2)
img = expand_image(a, 100)
return expand_image(b, 255, out=img)
img1 = expand_image(a, 255)
img2 = hitmiss_demo(a, [[0,0,0],[0,1,0],[1,1,1]], [[1,0,0],[0,0,0],[0,0,0]])
img3 = hitmiss_demo(a, [[0,0,0],[0,0,0],[1,1,1]], [[1,0,0],[0,1,0],[0,0,0]])
show_image(img1, img2, img3) 使用Hit和Miss運算的組合，可以實作複雜的圖像處理。例如文字識别中常用的細線化運算就可以用一系列的Hit和Miss運算實作。圖5顯示了細線化處理的效果，實作程式如下：
圖5 使用Hit和Miss進行細線化運算
def skeletonize(img):
h1 = np.array([[0, 0, 0],[0, 1, 0],[1, 1, 1]]) ❶
m1 = np.array([[1, 1, 1],[0, 0, 0],[0, 0, 0]])
h2 = np.array([[0, 0, 0],[1, 1, 0],[0, 1, 0]])
m2 = np.array([[0, 1, 1],[0, 0, 1],[0, 0, 0]])
hit_list = []
miss_list = []
for k in range(4): ❷
hit_list.append(np.rot90(h1, k))
hit_list.append(np.rot90(h2, k))
miss_list.append(np.rot90(m1, k))
miss_list.append(np.rot90(m2, k))
img = img.copy()
while True:
last = img
for hit, miss in zip(hit_list, miss_list):
hm = morphology.binary_hit_or_miss(img, hit, miss) ❸
# 從圖像中删除hit_or_miss所得到的白色點
img = np.logical_and(img, np.logical_not(hm)) ❹
# 如果處理之後的圖像和處理前的圖像相同，則結束處理
if np.all(img == last): ❺
break
return img
a = pl.imread("scipy_morphology_demo2.png")[:,:,0].astype(np.uint8)
b = skeletonize(a) ❶以圖6所示的兩個結構元素為基礎，構造4個形狀為(3, 3)的二維數組：h1、m1、h2、m2。其中h1和m1對應圖中左邊的結構元素，而h2和m2對應圖中右邊的結構元素，h1和h2對應白色結構元素，m1和m2對應黑色結構元素。❷将這些結構元素進行90°、180°、270°旋轉之後一共得到8個結構元素。
圖6 細線化算法的4個結構元素
❸依次使用這些結構元素進行Hit和Miss運算，❹并從圖像中删除運算所得到的白色像素，其效果就是依次從8個方向删除圖像的邊緣上的像素。❺重複運算直到沒有像素可删除為止。
二、圖像分割
下面以矩形區域識别為例，介紹如何使用measurements和morphology進行圖像區域分割。我們要抽取矩形資訊的圖像如圖7所示。這個圖像是二值圖像，其中矩形區域為白色，背景為黑色。但是由于它采用JPEG格式儲存，是以用pyplot.imread()讀取的是一個形狀為(高, 寬, 3)的三通道圖像。下面的程式使用其中的第0通道将其轉換成二值數組squares，将矩形區域設定為1，将背景設定為0。結果如圖7(左上)所示。
圖7 矩形區域分割算法各個步驟的輸出圖像
squares = pl.imread("suqares.jpg")
squares = (squares[:,:,0]＜200).astype(np.uint8) 由于許多矩形都有一些小凸起與鄰近的矩形連在一起，我們需要先将每塊矩形與其周圍的矩形分離出來。可以使用上節介紹的二值腐蝕函數morphology.binary_erosion()實作這一功能。不過這裡我們采用另外的方法。
morphology.distance_transform_cdt(image)計算二值圖像中每個像素到最近的黑色像素的距離，傳回一個儲存所有距離的數組。圖像上兩點之間的距離有很多定義方式，此函數預設采用切比雪夫距離。兩點之間的切比雪夫距離定義為其各坐标數值差的最大值，即DChess=max(｜x2-x1｜,｜y2-y1｜)。
下面調用distance_transform_cdt(squares)得到距離數組squares_dt，并繪制成圖。圖中顔色越紅的像素表示該點距離黑色背景越遠，而原圖中值為0的像素對應的距離為0，離黑色背景最遠的距離為27個像素。如果将距離數組當作圖像輸出，顯示效果如圖7(中上)所示。
from scipy.ndimage import morphology
squares_dt = morphology.distance_transform_cdt(squares)
print "各種距離值", np.unique(squares_dt)
各種距離值 [ 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] 隻需要選擇合适的門檻值将距離數組squares_dt轉換成二值數組，就可以實作矩形區域的分離，其效果和腐蝕算法類似。squares_core中每個矩形區域都縮得足夠小，以至于沒有任何兩塊區域之間有連通的路徑。其效果如圖7(右上)所示，可以看到所有的區域都和其相鄰的區域分割開了。
squares_core = (squares_dt＞8).astype(np.uint8) 下面調用measurements中的label()和center_of_mess()對各個獨立的白色區域進行染色，并計算各個區域的重心坐标。
labels()對二值圖像中每個獨立的白色區域使用唯一的整數進行填充，這相當于将每塊區域染上不同的“顔色”。這裡所謂的“顔色”是指為每個區域指定的一個整數，而不是指圖像中像素的真正顔色。
labels()傳回的結果數組可以用于計算各個區域的一些統計資訊。例如可以調用center_of_mass()計算每個區域的重心。其第一個參數為各個像素的權值(可以認為是每個像素的品質密度)，第二個參數為labels()輸出的染色數組，第三個參數為需要計算重心的區域的标簽清單。在下面的程式中，權值數組和染色數組相同，是以可以計算區域中所有白色像素的重心。
如果直接使用imshow()顯示labels()輸出的染色數組，将會得到一個區域顔色逐漸變化的圖像，這樣的圖像不利于肉眼識别各個區域。是以下面的程式用random_palette()為每個整數配置設定一個随機顔色。其結果如圖7(左下)所示，每個區域的重心使用白色小圓點表示。
from scipy.ndimage.measurements import label, center_of_mass
def random_palette(labels, count, seed=1):
np.random.seed(seed)
palette = np.random.rand(count+1, 3)
palette[0,:] = 0
return palette[labels]
labels, count = label(squares_core)
h, w = labels.shape
centers = np.array(center_of_mass(labels, labels, index=range(1, count+1)), np.int)
cores = random_palette(labels, count) 我們幾乎完成了區域識别的任務，但是每個矩形區域都比原始圖像中的要小一圈。下面将染色之後的矩形區域恢複到原始圖像中的大小。即對于原始圖像中為白色、腐蝕之後的圖像中為黑色的每個像素，将其顔色設定為染色數組中與之最近的區域的顔色。具體的計算步驟如下：
❶将腐蝕之後的圖像square_core反轉，并調用distance_transform_cdt()。這樣可以找到square_core中距離每個黑色像素最近的白色像素的坐标。為了讓其傳回坐标資訊，需要設定參數return_indices為True。由于這裡不需要距離資訊，是以可以設定參數return_distances為False。其傳回值index是一個形狀為(2, 高, 寬)的三維數組。index[0]為最近像素的第0軸(縱軸)坐标，index[1]為最近像素的第1軸(橫軸)坐标。
❷使用index[0]和index[1]擷取labels中對應坐标的顔色，得到near_labels。
❸建立一個布爾數組mask，其中的每個True值對應squares中為白色、squares_core中為黑色的像素。将labels複制一份到labels2，最後将mask中True對應的near_labels中的顔色值複制到labels2中。
圖7(中下)是使用random_palette()随機着色之後的結果。
index = morphology.distance_transform_cdt(1-squares_core,
return_distances=False,
return_indices=True) ❶
near_labels = labels[index[0], index[1]] ❷
mask = (squares - squares_core).astype(bool)
labels2 = labels.copy()
labels2[mask] = near_labels[mask] ❸
separated = random_palette(labels2, count)