第十章 數學形態學 10.1 導論 10.2 基礎概念 10.3 膨脹與侵蝕 10.4 開啟與關閉 10.5 型態與交離轉換 10.6 其他型態演算 10.7 MATLAB的bwmorph函數 10.8 灰階形態學 10.9 灰階形態學的應用

10.1 導論 數學形態學，簡稱形態學，專門用來處理、分析影像中的形狀。本章會介紹基礎形態學工具，並探討二元數位影像之範例，並延伸使至灰階影像的應用。

10.2 基礎概念 10.2.1 平移 設A為二元數位影像中一個像素集合，W＝（x，y）為某個向量座標。 Aω為集合A平移向量（x，y）結果。 定義:Aω＝{(a,b) + (x,y) : (a,b) ∈A}.

圖10.1平移向量 W＝（2 , 2）。A集合在x軸方向與y軸方向產生位移，位移量由w決定，x軸往下延伸，y軸往右延伸。 0 1 2 3 4 0 1 2 3 4 0 0 1 1 2 2 3 3 4 4 （圖10.1 平移） ● ●

10. 2. 1 平移反射標記為Â，以原點為中心反射： 定義:Â ＝ {(-x,-y) : (x,y) ∈ A)}. 如圖10 -3 -2 -1 0 1 2 3 -3 -2 -1 1 2 3 （圖10.2 反射） ○ ●

10.3 膨脹與侵蝕

10.3.1膨脹 像素集合A與B，使用B來”膨脹”A，記為A♁B 定義為：A♁B ＝∪ AX. X∈ B 代表每個點x ∈ B ，都使用這些座標對A平移。 膨脹也符合交換律 定義為 ：A♁B ＝B♁A 範例中 B(位移向量)＝{(0,0),(1,1),(-1,1),(1,-1),(-1,-1)}用來平移A之座標

如圖10.3，膨脹有放大效果。根據B座標不同，A ♁B會出現於不同位置。 1 2 3 4 5 1 2 3 4 -1 0 1 5 -1 6 0 7 1 A (0,0) B ● ●

下列為經過 B＝{(0,0),(1,1),(-1,1),(1,-1),(-1,-1)}平移A座標 藍色部份為A（0，0）A本身。 1 2 3 4 5 1 2 3 4 5 1 1 2 2 3 3 4 4 5 5 6 6 7 7 A（1，1） A（-1，1） ● ●

1 2 3 4 5 1 2 3 4 5 1 1 2 2 3 3 4 4 5 5 6 6 7 7 A（1，-1） A（-1，-1） ● ●

A ♁B可以將A的所有像素點(x,y)用一個B的複製替換，將B的(0,0)點置於(x,y)的位置。同樣也可以複製A。 1 2 3 4 5 1 2 3 4 5 6 7 (圖為A ♁ B) ●

根據B座標不同，A ♁B可能出現離A很遠的地方，如下圖 A位置相同但B位置不同，B座標為B={(7,3),(6,2),(6,4),(8,2)(8,4)} A ♁B=A(7,3) ∪A(6,2) ∪A(6,4) ∪A(8,2) ∪A(8,4). 1 2 3 4 1 2 3 4 5 6 7 8 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 10 11 12 13 14 ● ●

使用Matlab進行膨脹 膨脹範例 影像變得粗厚，正是膨脹作用。 (3*3矩陣進行膨脹) >> t=imread(‘text.tif’); >>sq=ones(3,3); >>td=imdilate(t,sq); >>subplot(1,2,1),imshow(t) >>subplot(1,2,2),imshow(td) 圖10.5(原圖) (膨脹後)

10.3.2侵蝕 定義:A Ө B 由像素點w=(x,y)構成，將B置於A之內，找到相對應的B原點(0,0)所在位置之像素。這些像素點變形成侵蝕。 侵蝕A Ө B是A的子集合。要看B原點的位置而 1 2 3 4 5 6 1 2 3 -1 0 1 4 -1 5 0 6 1 圖10.6 (A) (B) ● ●

下列圖10.6為侵蝕範例，左圖為原圖A，藍色部分為B，右圖為侵蝕結果，視B圓點位置決定。 1 2 3 4 5 6 1 2 3 4 5 6 1 1 2 2 3 3 4 4 5 5 6 6 ● ●

下列為不同的侵蝕點，侵蝕結果會包含在原始圖中。 1 2 3 4 5 6 1 2 3 4 5 6 1 1 2 2 3 3 4 4 5 5 6 6 ● ●

1 2 3 4 5 6 1 2 3 4 5 6 1 1 2 2 3 3 4 4 5 5 6 6 ● ●

下圖10.7是B不包含原點座標的侵蝕結果。空心圓為侵蝕結果。 侵蝕形狀與圖10.6相同但位置不同 0 1 2 3 4 0 1 2 3 4 5 6 0 -2 1 -1 2 0 3 1 4 2 5 3 6 4 B 5 6 A Ө B ○ ● ○ ●

使用Matlab進行侵蝕 侵蝕範例 如圖10.8影像變細，正是侵蝕預期的作用 (3*3十字矩陣進行侵蝕) sq=ones(3,3); c=imread(‘circbw.tif’); ce=imerode(c,sq); subplot(1,2,1),imshow(c) 圖10.8 (原圖) subplot(1,2,2),imshow(ce) (侵蝕後)

10.3.3應用:邊界偵測 若A為一影像，B為以原點為中心對稱分佈的小型結構元素，A邊界可以由下列幾種方法定義 (1) 內被邊界 A - (A Ө B) (2) 外部邊界(A ♁ B) – A (3)形態梯度 (A ♁B) - (A Ө B) 負號代表集合相減

下列範例圖10.9為(A), (B), (A♁B), (AӨB) 0 1 2 3 4 5 -1 0 1 2 3 4 5 6 0 (A) -1 1 0 2 1 3 2 4 3 5 4 0 1 2 3 4 5 5 0 (A Ө B) 6 1 2 -1 0 1 (A ♁ B) 3 -1 (B) 4 0 5 1 ● ● ● ●

下列為內部邊界,外部邊界與形態梯度。 『內部邊界』由A邊緣的像素組成。 『外部邊界』由A相鄰但位於A之外的像素組成。 『形態梯度』由內部邊界與外部邊界的總合。 -1 0 1 2 3 4 5 6 -1 0 1 2 3 4 5 6 -1 -1 0 0 1 1 2 2 3 3 4 4 5 5 6 6 (實點為內部邊界白點為外部邊界) (形態梯度) ○ ● ●

以下函數程式為稻穀影像與內部邊界 以圖10.10稻穀影像為例，設定閥值，求得二元 數位影像。 rice=imread(‘rice.tif’); r=rice>110; 再用下列指令可顯示原圖與內部邊界。 sq=ones(3,3); re=imerode(r,sq); (稻穀原圖) r_int=r&~re; subplot(1,2,1),imshow(r) subplot(1,2,2),imshow(r_int) (內部邊界)

以下指令為外部邊界與形態梯度 外部邊界與形態梯度也可以用類似方式求得。 re=imerode(r,sq); rice=imread(‘rice.tif’); r=rice>110; sq=ones(3,3); rd=imdilate(r,sq); 圖10.11 (外部邊界) r_ext=rd&~r; r_grad=rd&~re; subplot(1,2,1),imshow(r_ext) subplot(1,2,2),imshow(r_grad) (形態梯度)

10.4 開啟與關閉

10.4.1.開啟 就是先侵蝕後膨脹。如下圖10.12解釋 (B) 0 (A) (A Ө B) (A ○ B) ● A與結構元素B，用B來開啟A，符號為A ○ B。 -1 0 1 定義: A ○ B = (A Ө B) ♁ B -1 就是先侵蝕後膨脹。如下圖10.12解釋 (B) 0 1 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 (A) (A Ө B) (A ○ B) ● ● ● ●

10.4.2關閉 相對於開啟，關閉則是先膨脹再侵蝕，符號為A ● B 定義: A ● B = (A ♁ B) Θ B. 如下圖10.13解釋 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 1 2 3 4 5 6 0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 (A) (A ♁ B ) A ● B ● ● ●

開啟跟關閉可分別使用imopen與imclose執行。下面可以看到使用方形與十字結構元素的效果。

為imclose。

使用關閉，影像會完全黏合。使用斜角線結構元素黏合影像。 如圖10.14為關閉範例 t=imread('text.tif'); diag=[0 0 1; 0 1 0; 1 0 0] (斜角結構元素) Diag = 0 0 1 0 1 0 1 0 0 tc=imclose(t,diag); imshow(tc) (關閉範例) (關閉範例)

應用:去除雜訊 A ӨB會去除單一黑色像素，但會使空洞變大，連續膨脹兩次便可填充空洞((A ӨB) ♁B) ♁B. 第一次膨脹:恢復原始大小。 第二次膨脹:去除空洞。(但會放大影像中的物體) 恢復原始尺寸>>執行一次侵蝕(((A Ө B) ♁ B) ♁ B) Ө B. 內層兩個運算為開啟，外層為關閉。 所以去除雜訊方法為:先執行開啟後執行關閉((A ○ B) ● B). 這樣運算稱為:『型態濾波』。

加入10％散粒雜訊 如圖10.15(a)受到脈衝雜訊 干擾的二元雜訊影像。 c=imread(‘circles.tif’) ; x=rand(size(c)); d1=find(x<=0.05); d2=find(x>=0.95); c(d1)=0; c(d2)=1; imshow(c) (a) (二元雜訊影像)

濾波過程如下 如下圖10.15(b),(c) cf1=imclose(imopen(c,sq),sq); (使用方型結構元素) figure,imshow(cf1) cf2=imclose(imopen(c,cr),cr); (使用十字結構元素) figure,imshow(cf2) 結果方形結構元素有些鋸齒狀，若使用十字型結構元素則較平滑。 (b) (使用方形結構元素) (c)(使用十字結構元素)

10.5 型態交離轉換

假設要尋找3×3的方形，如下圖10.16(A)中央部分。 （1）影像A與B執行侵蝕A Ө B 型態交離轉換功能:影像中尋找特定形狀。 假設要尋找3×3的方形，如下圖10.16(A)中央部分。 （1）影像A與B執行侵蝕A Ө B (A) (B) ● ●

(2) 得結果如圖10.17，剩下兩個像素，因為影像A只有這兩個地方與B符合。 ●

(3) 如下圖10.18使用結構C侵蝕A的補數(C為3×3外部邊界圖) ● ● ●

使用如上一章節文字影像中”Corss-Correlation”中的連字部分符號。 產生結構元素程式如下: b1=ones(1,6); b2=[1 1 1 1 1 1 1 1;1 0 0 0 0 0 0 1;1 1 1 1 1 1 1 1]; tb1=erode(t,b1); tb2=erode(~t,b2); hit_or_miss=tb1&tb2; [x,y]=find(hit_or_miss==1)

10.6 其他形態演算法

10.6.1 區域填充 設一個區域為連通邊界畫分出來，如下圖10.21像素p為區域內的一個像素點。 ※要填充整個區域從p開始: (1) 使用十字結構元素B連續執行膨脹。 (2)執行下一次膨脹前與Ā做交集。 最後得到數列 {p} = X0,X1,X2….,Xk = Xk+1, 其中 Xn = (Xn-1 ♁ B) ∩Ā. ●

下圖10.22填充區域之過程 Xk∪A為填充區域，下圖顯示整個流程。 X0 = {p}, X1 = {p,1}, X2 = {p,1,2}… (Ā) ● 6 5 4 3 2 1 p

10.6.2連通部分 ● 如右圖10.23填充連通部分 連通部分也可以使用類似演算法。 4-連通部分使用十字結構元素。 8-連通部分使用方形結構元素。 像素p開始，產生一個集合的數列來填充。 X0 = {p},X1,X2,……, 其中 Xn = (Xn-1 ♁ B) ∩ A. 直到Xk = Xk-1，範例如右圖 十字結構元素無法穿越邊界。 (使用十字形) (使用方形) 2 1 p 5 4 3 2 1 p

進行MATLAB函數進行填充區域時，持續記錄兩個影像: (current與last)，last影像從區域的單一像素點p開始，current影像從膨脹開始。 Last  current, Current (current ♁B) ∩ Ā 給定B，執行: Imsilate(current,B)&~A. 可以用下列函數填充由邊界劃分的區域: sq=ones(3,3); n=imread('nicework.tif'); imshow(n),pixval on nb=n&~imerode(n,sq); figure,imshow(nb) nf=regfill(nb,[74,52],sq); figure,imshow(nf) 原圖(a) (b) 區域填充另一種表示法

以下函數可以來填充某個由邊界畫分的區域 (c) (d) (c)是區域填充結果 (d)是另一種表示方法 function out=regfill(im,pos,kernel); n=imread('nicework.tif'); current=zeros (size(im)); last=zeros(size(im)); last(pos(1),pos(2))=1; current=imdilate(last,kernel)&~im; while any(current(:)~=last(:)), last=current; end; out=current; (c) (d)

連通部分可用影像“nice work”來測試此函數，使用方形結構元素及一個較大的11*11結構元素。 (a) 為使用3*3方形 (b) 為使用11*11方形。 sq=ones(3,3); n=imread('nicework.tif'); sq2=ones(11,11); nc=components(n,[57,97],sq); imshow(nc) nc2=components(n,[57,97],sq2); figure,imshow(nc2) (a) 。 (b)

連通部份的函數和區域填充幾乎一模一樣，只要更改一行指令即可。 function out = components(im,pos,kernel) current=zeros (size(im)); last=zeros(size(im)); last(pos(1),pos(2))=1; current=imdilate(last,kernel)&im; while any(current(:)~=last(:)), last=current; end; out=current;

10.6.3 骨架化 二元物體的骨架式壓縮了物體大小與形狀後，所得到的直線與曲線的集合。 ※表10.1 建構骨架的運算 侵蝕 開啟 差值 A A ○ B A－(A ○ B) A Ө B (A Ө B) ○ B (A Ө B) －((A Ө B) ○ B) A Ө 2B (A Ө 2B) ○ B (A Ө 2B) －((A Ө 2B) ○ B) A Ө 3B (A Ө 3B) ○ B (A Ө 3B) －((A Ө 3B) ○ B) . A Ө kB (A Ө kB) ○ B (A Ө kB) －((A Өk B) ○ B)

範例如下圖10.28，骨架化，使用十字形結構元素。 (A) 原圖 ( A ○ B) 開啟 (A－(A ○ B))差值 ● ● ●

(A Ө B)侵蝕 ((A Ө B) ○ B)開啟 ((A Ө B) －((A Ө B) ○ B))差值 ● ● ●

(A Ө 2B)侵蝕 ((A Ө 2B) ○ B)開啟 ((A Ө 2B) －((A Ө 2B) ○ B))差值 ● ●

因為(A Ө 2B) ○ B已經等於零，便停止運算。骨架的第三欄中所有集合加起來就是最終骨架(圖10.29)。 ● 圖10.29(最終骨架) ● ● ●

若要用Matlab實現這種演算法，可寫成下列函數 如右下圖10.31是使用方型結構元素。 右方為函數指令 sq=ones(3,3); n=imread('nicework.tif'); nk=imskel_in(n,sq); imshow(nk) nk2=imskel_in(n,cr); figure,imshow(nk2) function skel = imskel (image,str) skel = zeros(size(image)); e = image; while (any(e(:))), o=imopen(e,str); nicework=skel|(e&~o); e=imerode(e,str); end

10.7 MATLAB的bwmorph函數 考慮一像素周圍之3*3鄰域。因為鄰域中像素只有兩種數值，所以2的9次方=512種可能的鄰域。 每種可能的鄰域狀態都對應一個從0到511的數值。 運用此方法，可以將膨脹定義如下: 一個數值為0的像素，若其鄰接的8個像素中，至少有一個白色的數值為1，則該像素變為1的像素。 若鄰接的8個像素中，至少有一個數值為0，則該像素變為0的像素。 優點:可以透過定義查詢表格來執行。如骨架(線條的粗細就是一個像素寬) 缺點:只能使用3*3鄰域。

10.8灰階形態學

二元侵蝕設B為3*3方形，數值為0 (A)為圖3.32的影像。 在影像A移動有下列步驟 1. 找出p的3*3鄰域Np。 2. 計算矩陣Np-B 3. 找出結果的最小值。 (A) (B) 1

由於B的數值皆為0，因此第二,第三步驟可以簡化為找出Np的最小值。 若鄰域中包含至少一個0，輸出為0。 若鄰域中包含至少一個1，則輸出為1。 例 0 0 0 0 1 1 0 0 1 1 若執行此運算，可得右圖 1

至於膨脹步驟和侵蝕類似 1. 找出p的3*3鄰域Np。 2. 計算矩陣Np+B 3.找出結果的最大值。 執行上列步驟可得右圖 1

若A為灰階影像，B為結構元素，則灰階影像可定義對影像中的各像素p執行下列步驟: 1. 將B置於A上方，(0,0) 點對準p。 2. 找出對應B形狀之p的鄰域Np。 3.找出數值min (Np-B)。 這些形成B點的集合稱”B點的定義域”，記做DB。 定義為:(A Ө B)(x,y) = min{A(x+s,y+t)-B(s,t),(s,t) ∈ DB}, (A♁B)(x,y) = max{A(x+s,y+t)+B(s,t),(s,t) ∈ DB},

灰階侵蝕與膨脹之間的關係 依照最大值與最小值的定義，若X與Y為兩矩陣 則:max{X+Y} = -min{-X-Y} 由於max{X+Y} 與(A♁B)相對應 min{X-Y} 與(A Ө B)相對應 則: A ♁ B = -(-A Ө B) A Ө B = -(A ♁ B) 或 -(A ♁ B) = -A Ө B -(A Ө B) = -A ♁ B

可以使用含數imerode與imdilate執行灰階侵蝕與膨脹。 首先須產生結構元素。 str=strel('arbitrary',ones(3,3),[1 2 3;4 5 6;7 8 9]) 在此設定strel的參數arbitrary。第一個矩陣設定鄰域，第二個矩陣提供數值。如此可進行測試 (測試)

對於膨脹，Matlab將結構元 素旋轉180度後執行。 必須旋轉str來獲得str2。 str2=strel('arbitrary',ones(3,3), [9 8 7;6 5 4;3 2 1]) imdilate(A,str2)

現在使用影像測試。 輸入下列函數: 可得右圖10.34(a),(b) c=imread('caribou.tif'); str=strel('square',5) cd=imdilate(c,str); ce=imerode(c,str); imshow(cd),figure,imshow(ce) a(膨脹) b(侵蝕)

開啟與關閉 如右圖10.35 灰階影像可以使用函數imopen與imclose，不過結構元素必須使用strel產生。 輸入下列函數: co=imopen(c,str); cc=imclose(c,str); imshow(co),figure,imshow(cc) (a)開啟 (b)關閉

10.9灰階型態學的應用

10.9.1 邊緣偵測 可以使用形態梯度 (A ♁B) - (A Ө B) 下面使用兩種結構元素來實驗 (3*3)與(5*5)如右圖10.36 str1=strel('square',3); str2=strel('square',5); ce1=imerode(c,str1); ce2=imerode(c,str2); cd1=imdilate(c,str1); cd2=imdilate(c,str2); cg1=imsubtract(cd1,ce1); cg2=imsubtract(cd2,ce2); imshow(cg1),figure,imshow(cg2) (a) 3*3結構元素 (b) 5*5結構元素

10.9.2 去除雜訊 使用形態濾波來去除雜訊: 先開啟後關閉 ((A ○ B) ● B). 如右圖10.37 型態濾波對去除高斯雜訊效果不太好。 cn=imnoise(c,'salt & pepper'); cf=imclose(imopen(cn,str1),str); imshow(cn),figure,imshow(cf) (a) (b)

報告結束 謝謝各位~