腐蝕膨脹-算法

本篇文章要分享的是基於MATLAB的腐蝕膨脹算法實現，腐蝕膨脹是形態學圖像處理的基礎，腐蝕在二值圖像的基礎上做“收縮”或“細化”操作，膨脹在二值圖像的基礎上做“加長”或“變粗”的操作。

什么是二值圖像呢？把一幅圖片看做成一個二維的數組，那么二值圖像是一個只有0和1的邏輯數組，我們前面Sobel邊緣檢測后的圖像輸出邊緣效果，設置個閾值，大於閾值輸出為1，小於閾值輸出為0，最后輸出就是一幅二值圖像了。

腐蝕

腐蝕是一種消除邊界點，使邊界向內部收縮的過程。可以用來消除小且無意義的物體。用3X3的結構元素，掃描圖像的每一個像素，用結構元素與其覆蓋的二值圖像做“與”操作，如果都為1，結果圖像的該像素為1。否則為0。結果會使二值圖像小一圈。

有一個形象的比喻來可以說明該運算，用0表示蛀蟲，1表示大米。蛀蟲腐蝕大米的過程便是

 

 

 

腐蝕

如圖所示，對於一個像素矩陣而言，只要有蛀蟲（0）的存在，大米（1）就會被腐蝕掉了，即使只存在一個蛀蟲（0），但是還是會被蛀蟲腐蝕完畢，最后一幅圖上面由於沒有蛀蟲（0）所以大米完好無損。

關於算法的實現，可以用下式子來表示，即3x3像素的運算：

P = P11 & P12 & P13 & P21 & P22 & P23 & P31 & P32 & P33

在FPGA中，為了通過面積去換速度，我們將上式改變如下：　　　　　　　　　　　　　　　　　　　　　　　

P1 = P11 & P12 & P13

P2 = P21 & P22 & P23

P3 = P31 & P32 & P33

P = P1 & P2 & P3

MATLAB中可以直接寫一個按位或運算。

膨脹

膨脹是將與物體接觸的所有背景點合並到該物體中，使邊界向外部擴張的過程。可以用來填補物體中的空洞。用3X3的結構元素，掃描圖像的每一個像素，用結構元素與其覆蓋的二值圖像做“與”操作，如果都為0，結果圖像的該像素為0,。否則為1。結果使二值圖像擴大一圈。

先腐蝕后膨脹的過程稱為開運算。用來消除小物體、在纖細點處分離物體、平滑較大物體的邊界的同時並不明顯的改變其面積。先膨脹后腐蝕的過程稱為比運算，用來填充物體內細小空間、連接鄰近物體、平滑其邊界的同時並不明顯改變其面積。

膨脹算法用最簡單的比喻來描述：0表示害蟲，1表示青蛙，青蛙吃了害蟲表示膨脹運算，我們用3*3像素陣列來解釋：

 

 

 

膨脹

如圖所示，圖左只有害蟲（0），所以害蟲都活着，中間那個圖，雖然只有一個害蟲，但是還是會被青蛙全部吃掉，最右邊的那幅圖，都是青蛙，所以青蛙始終是青蛙。

關於算法的實現，可以用下式子來表示，即3x3像素的運算：

P = P11 | P12 | P13 | P21 | P22 | P23 | P31 | P32 | P33

在HDL中，為了通過面積去換速度，我們將上式改變如下：　　　　　　　　　　　　　　　　　　　　　　　　　　　　　

P1 = P11 | P12 | P13

P2 = P21 | P22 | P23

P3 = P31 | P32 | P33

P = P1 | P2 | P3

MATLAB中可以直接寫一個按位與運算。

開運算閉運算

先腐蝕后膨脹叫開運算，開運算的作用是清除圖像邊緣周圍非邊緣的細小的點。先膨脹后腐蝕為閉運算，閉運算的作用是清除圖像內部的空洞，

如果我們的目標物體外面有很多無關的小區域，就用開運算去除掉；如果物體內部有很多小黑洞，就用閉運算填充掉。

MATLAB邏輯運算函數

bitand(), 對十進制數進行逐位邏輯與運算：先將十進制數轉換成二進制數，然后逐位與運算，其運算結果轉換為十進制。

bitor(), 對十進制數進行逐位邏輯或運算：先將十進制數轉換成二進制數，然后逐位與運算，其運算結果轉換為十進制。

MATLAB代碼實現

%%%imopen Erosion_Dilation 
clc,clear,close all,clear all

img = imread('1.png');
img_gray = rgb2gray(img);
 figure;
subplot(2,2,1),imshow(img),title('灰度圖像') 
subplot(2,2,2),imshow(img_gray),title('灰度圖像')
hold on;
%%中值濾波
 %Median Filter 
 imgn = imnoise(img_gray,'salt & pepper',0.02); 
 subplot(2,2,3),imshow(imgn),title('椒鹽噪聲圖像'); 

 Median_Img = img_gray; 
 [ROW,COL] = size(Median_Img);
 for r = 2:ROW-1 
     for c = 2:COL-1 
         median3x3 =[imgn(r-1,c-1)    imgn(r-1,c) imgn(r-1,c+1) 
                     imgn(r,c-1)      imgn(r,c)      imgn(r,c+1) 
                     imgn(r+1,c-1)      imgn(r+1,c) imgn(r+1,c+1)]; 
         sort1 = sort(median3x3, 2, 'descend');
         sort2 = sort([sort1(1), sort1(4), sort1(7)], 'descend'); 
         sort3 = sort([sort1(2), sort1(5), sort1(8)], 'descend'); 
         sort4 = sort([sort1(3), sort1(6), sort1(9)], 'descend'); 
         mid_num = sort([sort2(3), sort3(2), sort4(1)], 'descend'); 
         Median_Img(r,c) = mid_num(2); 
     end 
 end 
  
subplot(2,2,4),imshow(Median_Img),title('中值濾波圖像')

%%Sobel算子邊緣檢測
 Median_Img = double(Median_Img);
 Sobel_Img = zeros(ROW,COL);
 Sobel_Threshold = 80;
for r = 2:ROW-1
     for c = 2:COL-1
        Sobel_x = Median_Img(r-1,c+1) + 2*Median_Img(r,c+1) + Median_Img(r+1,c+1) - Median_Img(r-1,c-1) - 2*Median_Img(r,c-1) - Median_Img(r+1,c-1);
        Sobel_y = Median_Img(r-1,c-1) + 2*Median_Img(r-1,c) + Median_Img(r-1,c+1) - Median_Img(r+1,c-1) - 2*Median_Img(r+1,c) - Median_Img(r+1,c+1);
        %Sobel_Num = abs(Sobel_x) + abs(Sobel_y);
         Sobel_Num = sqrt(Sobel_x^2 + Sobel_y^2);
        if(Sobel_Num > Sobel_Threshold)
             Sobel_Img(r,c)=1;
         else
             Sobel_Img(r,c)=0;
         end
     end
end
figure
subplot(2,2,1),imshow(Sobel_Img),title('Sobel算子圖像')
BW1=edge(img_gray,'sobel'); %用Sobel算子進行邊緣檢測
subplot(2,2,2),imshow(BW1),title('Sobel算子圖像')

%%腐蝕算法
Erosion_img = zeros(ROW,COL); 
for r = 2:ROW-1
     for c = 2:COL-1
        %and1 = bitand(Sobel_Img(r-1,c-1) + bitand(Sobel_Img(r-1,c) + Sobel_Img(r-1,c+1)));
        %and2 = bitand(Sobel_Img(r  ,c-1) + bitand(Sobel_Img(r  ,c) + Sobel_Img(r  ,c+1)));
        %and3 = bitand(Sobel_Img(r+1,c-1) + bitand(Sobel_Img(r+1,c) + Sobel_Img(r+1,c+1)));
        %Erosion_img(r,c) = bitand(and1,bitand(and2,and3));
        and1 = Sobel_Img(r-1,c-1) & Sobel_Img(r-1,c) & Sobel_Img(r-1,c+1);
        and2 = Sobel_Img(r  ,c-1) & Sobel_Img(r  ,c) & Sobel_Img(r  ,c+1);
        and3 = Sobel_Img(r+1,c-1) & Sobel_Img(r+1,c) & Sobel_Img(r+1,c+1);
        Erosion_img(r,c) = and1 & and2 & and3;
     end
end
subplot(2,2,3),imshow(Erosion_img),title('腐蝕圖像')
%%膨脹算法
Dilation_img = zeros(ROW,COL); 
for r = 2:ROW-1
     for c = 2:COL-1
        %and1 = bitand(Sobel_Img(r-1,c-1) + bitand(Sobel_Img(r-1,c) + Sobel_Img(r-1,c+1)));
        %and2 = bitand(Sobel_Img(r  ,c-1) + bitand(Sobel_Img(r  ,c) + Sobel_Img(r  ,c+1)));
        %and3 = bitand(Sobel_Img(r+1,c-1) + bitand(Sobel_Img(r+1,c) + Sobel_Img(r+1,c+1)));
        %Erosion_img(r,c) = bitand(and1,bitand(and2,and3));
        and1 = Sobel_Img(r-1,c-1) | Sobel_Img(r-1,c) | Sobel_Img(r-1,c+1);
        and2 = Sobel_Img(r  ,c-1) | Sobel_Img(r  ,c) | Sobel_Img(r  ,c+1);
        and3 = Sobel_Img(r+1,c-1) | Sobel_Img(r+1,c) | Sobel_Img(r+1,c+1);
        Dilation_img(r,c) = and1 | and2 | and3;
     end
end
subplot(2,2,4),imshow(Dilation_img),title('膨脹圖像')

%%開運算-先腐蝕后膨脹算法
Erosion_Dilation_img = zeros(ROW,COL); 
for r = 2:ROW-1
     for c = 2:COL-1
        %and1 = bitand(Sobel_Img(r-1,c-1) + bitand(Sobel_Img(r-1,c) + Sobel_Img(r-1,c+1)));
        %and2 = bitand(Sobel_Img(r  ,c-1) + bitand(Sobel_Img(r  ,c) + Sobel_Img(r  ,c+1)));
        %and3 = bitand(Sobel_Img(r+1,c-1) + bitand(Sobel_Img(r+1,c) + Sobel_Img(r+1,c+1)));
        %Erosion_img(r,c) = bitand(and1,bitand(and2,and3));
        and1 = Erosion_img(r-1,c-1) | Erosion_img(r-1,c) | Erosion_img(r-1,c+1);
        and2 = Erosion_img(r  ,c-1) | Erosion_img(r  ,c) | Erosion_img(r  ,c+1);
        and3 = Erosion_img(r+1,c-1) | Erosion_img(r+1,c) | Erosion_img(r+1,c+1);
        Erosion_Dilation_img(r,c) = and1 | and2 | and3;
     end
end
figure
subplot(1,2,1),imshow(Erosion_Dilation_img),title('先腐蝕在膨脹圖像')
%%閉運算-先膨脹后腐蝕算法
Dilation_Erosion_img = zeros(ROW,COL); 
for r = 2:ROW-1
     for c = 2:COL-1
        %and1 = bitand(Sobel_Img(r-1,c-1) + bitand(Sobel_Img(r-1,c) + Sobel_Img(r-1,c+1)));
        %and2 = bitand(Sobel_Img(r  ,c-1) + bitand(Sobel_Img(r  ,c) + Sobel_Img(r  ,c+1)));
        %and3 = bitand(Sobel_Img(r+1,c-1) + bitand(Sobel_Img(r+1,c) + Sobel_Img(r+1,c+1)));
        %Erosion_img(r,c) = bitand(and1,bitand(and2,and3));
        and1 = Dilation_img(r-1,c-1) & Dilation_img(r-1,c) & Dilation_img(r-1,c+1);
        and2 = Dilation_img(r  ,c-1) & Dilation_img(r  ,c) & Dilation_img(r  ,c+1);
        and3 = Dilation_img(r+1,c-1) & Dilation_img(r+1,c) & Dilation_img(r+1,c+1);
        Dilation_Erosion_img(r,c) = and1 & and2 & and3;
     end
end
subplot(1,2,2),imshow(Dilation_Erosion_img),title('先膨脹在腐蝕圖像')