1. 程式人生 > >影象處理(2)空間濾波

影象處理(2)空間濾波

  1. 濾波原理
    濾波是為了消除噪聲影響,提高影象質量,濾波方法的選擇取決於噪聲的特性。濾波包括空間濾波和頻率域濾波,在數學上二者可以相互轉化。空間濾波本質上是通過臨近畫素的微操作來實現噪聲降低,比如線性平滑濾波,就是通過臨近畫素的加權或者不加權平均來消減噪聲強度,這對於方差較小較穩定的噪聲效果較好。用公式來表示空間濾波就是:
    這裡寫圖片描述
    實際上就是通過一個w模板在影象上移動,然後進行相關性運算。也可以進行卷積運算,但是兩者在處理方法上沒有區別,卷積運算是將模板在影象水平和豎直方向進行翻轉後再進行相關操作。
  2. 線性濾波
    線性濾波有均值濾波,也有加權濾波。模板如下:
    這裡寫圖片描述
    MATLAB程式碼:
function filtered_image=linear_filtering(f,w,filtering_mode,boundary_options,size_options)

[row_f,col_f]=size(f);
[row_w,col_w]=size(w);

%expand the f iamge
f_exp=zeros(row_f+row_w-1,col_f+col_w-1);


f_sum=0;

%replicate the outer image data
if strcmp(boundary_options,'repllicate')
    for i
=1:row_f+row_w-1 for j=1:col_f+col_w-1 if (i<row_w/2) && (j<col_f+col_w/2 && j>col_w/2) f_exp(i,j)=f(1,j-floor(col_w/2)); elseif (i>row_f+row_w/2) && (j<col_f+col_w/2 && j>col_w/2) f_exp(i,j)=f(row_f,j
-floor(col_w/2)); elseif (j<col_w/2) && (i<row_f+row_w/2 && i>row_w/2) f_exp(i,j)=f(i-floor(row_w/2),j); elseif (j>col_f+col_w/2) && (i<row_f+row_w/2 && i>row_w/2) f_exp(i,j)=f(i-floor(row_w/2),col_f); elseif i<row_f+row_w/2 && i>row_w/2 && j<col_f+col_w/2 && j>col_w/2 f_exp(i,j)=f(i-floor(row_w/2),j-floor(col_w/2)); end end end %mirror of the image data in the outer space elseif strcmp(boundary_options,'symmetric') for i=1:row_f+row_w-1 for j=1:col_f+col_w-1 if (i<row_w/2) && (j<col_f+col_w/2 && j>col_w/2) f_exp(i,j)=f(row_w-1-i,j-floor(col_w/2)); elseif (i>row_f+row_w/2) && (j<col_f+col_w/2 && j>col_w/2) f_exp(i,j)=f(2*row_f+floor(row_w/2)-i,j-floor(col_w/2)); elseif (j<col_w/2) && (i<row_f+row_w/2 && i>row_w/2) f_exp(i,j)=f(i-floor(row_w/2),col_w-1-j); elseif (j>col_f+col_w/2) && (i<row_f+row_w/2 && i>row_w/2) f_exp(i,j)=f(i-floor(row_w/2),2*col_f+floor(col_w/2)-j); elseif i<row_f+row_w/2 && i>row_w/2 && j<col_f+col_w/2 && j>col_w/2 f_exp(i,j)=f(i-floor(row_w/2),j-floor(col_w/2)); end end end %extend the image as 2D periodic function elseif strcmp(boundary_options,'circular') for i=1:row_f+row_w-1 for j=1:col_f+col_w-1 if (i<row_w/2) && (j<col_f+col_w/2 && j>col_w/2) f_exp(i,j)=f(row_f-i,j-floor(col_w/2)); elseif (i<row_w/2) && (j<col_w/2) f_exp(i,j)=f(row_f-i,col_f-j); elseif (i>row_f+row_w/2) && (j<col_f+col_w/2 && j>col_w/2) f_exp(i,j)=f(i-row_f,j-floor(col_w/2)); elseif (i>row_f+row_w/2) && (j<col_w/2) f_exp(i,j)=f(i-row_f,col_f-j); elseif (i>row_f+row_w/2) && (j>col_f+col_w/2) f_exp(i,j)=f(i-row_f,j-col_f); elseif (i<row_w/2) && (j>col_f+col_w/2) f_exp(i,j)=f(row_f-i,j-col_f); elseif (j<col_w/2) && (i<row_f+row_w/2 && i>row_w/2) f_exp(i,j)=f(i-floor(row_w/2),col_f-j); elseif (j>col_f+col_w/2) && (i<row_f+row_w/2 && i>row_w/2) f_exp(i,j)=f(i-floor(row_w/2),j-col_f); elseif i<row_f+row_w/2 && i>row_w/2 && j<col_f+col_w/2 && j>col_w/2 f_exp(i,j)=f(i-floor(row_w/2),j-floor(col_w/2)); end end end end filtered_expand_image=f_exp; %filtering is done by correlation if strcmp(filtering_mode,'corr') for i=floor(row_w/2):floor(row_f+row_w/2)-1 for j=floor(col_w/2):floor(col_f+col_w/2)-1 for k=-floor(row_w/2):1:floor(row_w/2) for l=-floor(col_w/2):1:floor(col_w/2) f_sum=f_sum+w(k+floor(row_w/2)+1,l+floor(col_w/2)+1)*f_exp(i+k+1,j+l+1); end end filtered_expand_image(i,j)=f_sum; f_sum=0; end end %filtering is done by convolution elseif strcmp(filtering_mode,'conv') for i=floor(row_w/2):floor(row_f+row_w/2)-1 for j=floor(col_w/2):floor(col_f+col_w/2)-1 for k=-floor(row_w/2):1:floor(row_w/2) for l=-floor(col_w/2):1:floor(col_w/2) f_sum=f_sum+w(k+floor(row_w/2)+1,l+floor(col_w/2)+1)*f_exp(i-k+1,j-l+1); end end filtered_expand_image(i,j)=f_sum; f_sum=0; end end end %the output image is the same as the extanded image if strcmp(size_options,'full') filtered_image=filtered_expand_image; %the same as the input image elseif strcmp(size_options,'same') filtered_image=filtered_expand_image(floor(row_w/2)+1:floor(row_f+row_w/2),... floor(col_w/2)+1:floor(col_f+col_w/2)); end

頂層函式為呼叫處理函式並顯示影象:

function filtering_top


image=imread('./image/B.jpg');

%transformed to gray 
[d1,d2,d3]=size(image);
if d3==3
    R=image(:,:,1);
    G=image(:,:,2);
    B=image(:,:,3);
    gray=round((1/3)*(R+G+B));
elseif d3==2
    R=image(:,:,1);
    G=image(:,:,2);
    gray=round((1/2)*(R+G));
else
    gray=image;
end

figure('color','w');
imshow(gray,[0,255]);

%gray=double(gray)+20*randn(size(gray));
gray=imnoise(gray,'salt & pepper');

figure('color','w');
imshow(gray,[0,255]);

%linear filtering
%3X3
%w=(1/9)*[1,1,1;1,1,1;1,1,1];

%5X5
w=(1/25)*ones(5,5);

%filtering
filtered_gray=linear_filtering(gray,w,'corr','circular','same');

%w=[5,5];
%filtered_gray=stistical_filtering(gray,w,'circular','same');



figure('color','w');
imshow(filtered_gray,[0,255]);

產生20%的加性噪聲和椒鹽噪聲,線性濾波器可以起到平滑的效果,模板越大平滑效果越明顯,但是同時也模糊了影象。效果如圖:

這裡寫圖片描述
圖2 (a)原始影象 (b)20%加性噪聲 (c)平滑均值濾波3x3 (d) 平滑均值濾波5x5 (e)椒鹽噪聲 (f)對椒鹽噪聲用平滑均值濾波器3x3效果

  1. 統計排序濾波
    統計排序即是在一個區域對畫素排序,選擇其中某個順序的畫素作為值。比如中值濾波,是選擇中間值。MATLAB程式包含stistical_filtering.m和rank.m,後者是對模板內資料進行排序。
function filtered_image=stistical_filtering(f,w,boundary_options,size_options)

[row_f,col_f]=size(f);

%[row_w,col_w]=w, the delt region
row_w=w(1);
col_w=w(2);

%expand the f iamge
f_exp=zeros(row_f+row_w-1,col_f+col_w-1);

%replicate the outer image data
if strcmp(boundary_options,'repllicate')
    for i=1:row_f+row_w-1
        for j=1:col_f+col_w-1
            if (i<row_w/2) && (j<col_f+col_w/2 && j>col_w/2) 
                f_exp(i,j)=f(1,j-floor(col_w/2));

            elseif (i>row_f+row_w/2) && (j<col_f+col_w/2 && j>col_w/2)
                f_exp(i,j)=f(row_f,j-floor(col_w/2));

            elseif (j<col_w/2) && (i<row_f+row_w/2 && i>row_w/2)
                f_exp(i,j)=f(i-floor(row_w/2),j);

            elseif (j>col_f+col_w/2) && (i<row_f+row_w/2 && i>row_w/2)
                f_exp(i,j)=f(i-floor(row_w/2),col_f);

            elseif i<row_f+row_w/2 && i>row_w/2 && j<col_f+col_w/2 && j>col_w/2
                f_exp(i,j)=f(i-floor(row_w/2),j-floor(col_w/2));

            end


        end
    end

%mirror of the image data in the outer space    
elseif strcmp(boundary_options,'symmetric')
    for i=1:row_f+row_w-1
        for j=1:col_f+col_w-1
            if (i<row_w/2) && (j<col_f+col_w/2 && j>col_w/2) 
                f_exp(i,j)=f(row_w-1-i,j-floor(col_w/2));

            elseif (i>row_f+row_w/2) && (j<col_f+col_w/2 && j>col_w/2)
                f_exp(i,j)=f(2*row_f+floor(row_w/2)-i,j-floor(col_w/2));

            elseif (j<col_w/2) && (i<row_f+row_w/2 && i>row_w/2)
                f_exp(i,j)=f(i-floor(row_w/2),col_w-1-j);

            elseif (j>col_f+col_w/2) && (i<row_f+row_w/2 && i>row_w/2)
                f_exp(i,j)=f(i-floor(row_w/2),2*col_f+floor(col_w/2)-j);

            elseif i<row_f+row_w/2 && i>row_w/2 && j<col_f+col_w/2 && j>col_w/2
                f_exp(i,j)=f(i-floor(row_w/2),j-floor(col_w/2));

            end


        end
    end


%extend the image as 2D periodic function
elseif strcmp(boundary_options,'circular')
    for i=1:row_f+row_w-1
        for j=1:col_f+col_w-1
            if (i<row_w/2) && (j<col_f+col_w/2 && j>col_w/2) 
                f_exp(i,j)=f(row_f-i,j-floor(col_w/2));

            elseif (i<row_w/2) && (j<col_w/2)
                f_exp(i,j)=f(row_f-i,col_f-j);

            elseif (i>row_f+row_w/2) && (j<col_f+col_w/2 && j>col_w/2)
                f_exp(i,j)=f(i-row_f,j-floor(col_w/2));

            elseif (i>row_f+row_w/2) && (j<col_w/2)
                f_exp(i,j)=f(i-row_f,col_f-j);

            elseif (i>row_f+row_w/2) && (j>col_f+col_w/2)
                f_exp(i,j)=f(i-row_f,j-col_f);

            elseif (i<row_w/2) && (j>col_f+col_w/2)
                f_exp(i,j)=f(row_f-i,j-col_f);

            elseif (j<col_w/2) && (i<row_f+row_w/2 && i>row_w/2)
                f_exp(i,j)=f(i-floor(row_w/2),col_f-j);

            elseif (j>col_f+col_w/2) && (i<row_f+row_w/2 && i>row_w/2)
                f_exp(i,j)=f(i-floor(row_w/2),j-col_f);

            elseif i<row_f+row_w/2 && i>row_w/2 && j<col_f+col_w/2 && j>col_w/2
                f_exp(i,j)=f(i-floor(row_w/2),j-floor(col_w/2));

            end


        end
    end


end


filtered_expand_image=f_exp;
rank_vector=zeros(1,row_w*col_w);


%rank the data and get the middle data as the image pixel
for i=floor(row_w/2)+1:row_f+floor(row_w/2)
    for j=floor(col_w/2)+1:col_f+floor(col_w/2)
       vector_n=0;
       for k=-floor(row_w/2):1:floor(row_w/2)
           for l=-floor(col_w/2):1:floor(col_w/2)

               vector_n=vector_n+1;
               rank_vector(vector_n)=f_exp(i+k,j+l);

           end
       end

       middle_data=rank(rank_vector);
       filtered_expand_image(i,j)=middle_data;

    end  
end



%the output image is the same as the extanded image
if strcmp(size_options,'full')

    filtered_image=filtered_expand_image;
%the same as the input image
elseif strcmp(size_options,'same')

    filtered_image=filtered_expand_image(floor(row_w/2)+1:floor(row_f+row_w/2),...
        floor(col_w/2)+1:floor(col_f+col_w/2));
end
function middle_data=rank(vector)

row=length(vector);

for i=1:row
    for k=2:row
        if(vector(i)>vector(k))
            vector_temp=vector(i);
            vector(i)=vector(k);
            vector(k)=vector_temp;
        end

    end
end

middle_data=vector(floor(row/2)+1);

統計排序濾波比線性濾波對椒鹽噪聲的濾波效果好,但是對加性隨機噪聲濾波效果一般。

這裡寫圖片描述
圖3 (a)使用中值濾波3x3對加性噪聲濾波 (b)5x5中值濾波 (c)5x5中值濾波對椒鹽噪聲

參考文獻:
[1]: 數字影象處理,Rafael. C. Gonzalez, Richard. E. woods.
[2]: image processing using matlab, Rafael. C. Gonzalez.