[数字图像处理]The Simulation and Application of Frequency Domain Filters（自己实现常用的滤波器）

最新推荐文章于 2022-05-19 09:15:44 发布

ourhonor

最新推荐文章于 2022-05-19 09:15:44 发布

阅读量750

点赞数 1

分类专栏：数字图像处理文章标签： matlab

本文链接：https://blog.csdn.net/zzx2015/article/details/89498297

版权

数字图像处理专栏收录该内容

4 篇文章 0 订阅

订阅专栏

Introduction:

In this lab, I review the principles and properties of three kinds of filters (Sobel filter, Gaussian low pass and high pass, Butterworth notch filters). In order to make a better understand of these three filters, I apply these methods into different figures and compare the results of different filters.

Case1: Sobel filter

Principle of Sobel filter:

Sobel operator is mainly used for edge detection. It is a discrete difference operator used to calculate the approximate value of image brightness function. The core of Sobel edge detection lies in the convolution of pixel matrix. Convolution is very important for digital image processing, and many image processing algorithms are realized by convolution. The essence of convolution operation is the process of weighted summation of pixel values in the designated image region. The calculation process is that each pixel value in the image region is multiplied by each element of the convolution template respectively, and the sum of the convolution results is the result of convolution operation.

The convolution formula of the matrix is as follows:

The convolution operation between window M of 3x3 and convolution template C is as follows:

Gx and Gy are convolution factors of sobel, and the two factors are convolved with the original image as follows:

The notation A represents the original image:

Get the horizontal and longitudinal gray value Gx and Gy of each point in the image. Finally, the following formula is used to calculate the size of gray scale change:

But usually, for efficiency, you don't use the approximation of the square root, although we lose precision:

Pseudocode of algorithm:

calculate the product of Gx, Gy and each row of the template.
the convolution of two 3x3 matrices is the multiplication and addition of each row and each column.
Gx and Gy after the calculation of 3*3 template.
the Gx ^ 2 + Gy ^ 2 square root or directly to the Gx and Gy sum after taking absolute value.
set a threshold value, and the pixel value after operation is greater than the threshold value, and the output is all 1, and less than the threshold value, and the output is all 0.

Matlab codes:

clc;

clear all;

close all;

Image = imread('Q3_1.tif');%

Sobel2_11711118(Image);

function Sobel2_11711118(inputImage)

[ROW,COL, DIM] = size(inputImage);

Gray_data = inputImage;

size1=size(inputImage);

image=zeros(2+size1(1),2+size1(2));

for a=1:size1(1)

for b=1:size1(2)

image(a,b)=((-1)^(a+b))*inputImage(a,b);

end

F=fft2(image);

figure,imshow(log(abs(F)),[ ]);

%% Median Filter

imgn=Gray_data;

Median_Img = Gray_data;

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

% figure;

% imshow(Median_Img);

%% Sobel_Edge_Detect

Median_Img = double(Median_Img);%Convert image data to double to avoid loss of precision

Sobel_Threshold = 150;

Sobel_Img = zeros(ROW,COL);

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)=0;

else

Sobel_Img(r,c)=255;

end

size2=size(Sobel_Img);

image2=zeros(2+size2(1),2+size2(2));

for a=1:size1(1)

for b=1:size1(2)

image2(a,b)=((-1)^(a+b))*Sobel_Img(a,b);

end

F2=fft2(image2);

figure

imshow(log(abs(F2)),[ ]);

title('output')

figure;

imshow(Sobel_Img);

end

Result and analysis:

Because the origin result is terrible and not clearly, we add the step called Median Filter and we can find that the result become much better than before.

This figure is the frequency domain of input figure.

This figure shows that the frequency domain of output figure. We can find that the noise is much less than before.

this figure is we finally get. We can see this figure is much better than the result we only do Sobel filter.

In another way, we can also do the Sobel filter:

The result and codes are in the following:

clc;

clear all;

close all;

imputImg = imread('Q3_1.tif');

Sobel1_11711118(imputImg);

function Sobel1_11711118(img);

size1=size(img);

image=zeros(2+size1(1),2+size1(2));

for a=1:size1(1)

for b=1:size1(2)

image(a,b)=((-1)^(a+b))*img(a,b);

end

F=fft2(image);

figure,imshow(log(abs(F)),[ ]);

title('input')

%blur_img = imgaussfilt(img);

blur_img=img;

% masks

gx = [-1/2, 0, 1/2; -1, 0, 1; -1/2, 0, 1/2];

gy = gx';

% convolute

% same means the size of the new image is the same; By adding zero paddings

% before calculating and deleting them afterwards

conv_img_x = conv2(blur_img, gx, 'same');

conv_img_y = conv2(blur_img, gy, 'same');

conv_img = abs(conv_img_x) + abs(conv_img_y);

% normalize

norm_img = normalize(conv_img);

figure

imshow(norm_img )

title('output')

end

Case2: Gaussian low pass and high pass

Principle of Sobel filter:

Gaussian filtering is the process of weighted average of the whole image. The value of each pixel is obtained by weighted average of itself and other pixel values in the neighborhood. The specific operation of gaussian filtering is to scan every pixel in the image with a template (or convolution, mask) and replace the value of the center pixel of the template with the weighted average gray value of the pixel in the neighborhood determined by the template.

The expression of Gaussian function is:

For a 3*3 template:

(x, y) is the template coordinate, then w1=h (-1,1), w2=h (-1,0) ...Generate the corresponding 3*3 filter template; Sigma is the standard deviation.

Pseudocode of algorithm:

Read the image
Generate the kernel of Gaussian
Convolution
Get the new picture

Matlab codes:

clc;close all;

GaussianLow_11711118(30);

GaussianLow_11711118(60);

GaussianLow_11711118(160);

GaussianHigh_11711118(30);

GaussianHigh_11711118(60);

GaussianHigh_11711118(160);

function GaussianLow_11711118(setD0)

a=setD0;

I=imread('Q3_2.tif');

I=im2double(I);

M=2*size(I,1);

N=2*size(I,2); %the size of filter

u=-M/2:(M/2-1);

v=-N/2:(N/2-1);

[U,V]=meshgrid(u,v);

D=sqrt(U.^2+V.^2);

D0=setD0;

H=exp(-(D.^2)./(2*(D0^2))); %generate the Guassian filter

J=fftshift(fft2(I,size(H,1),size(H,2)));

G=J.*H;

L=ifft2(fftshift(G));

L=L(1:size(I,1),1:size(I,2));

figure;

imshow(L);

title(['Result of GLPF with D0 = ',num2str(a)]);

end

function GaussianHigh_11711118(setD0)

a=setD0;

img=imread('Q3_2.tif');

% img=rgb2gray(img);

% figure,imshow(img);

g=fft2(double(img));

%%fourier transform

g=fftshift(g);

%centralization

%imshow(abs(g),[])

[N1,N2]=size(g);

n=2;

d0=setD0;

n1=fix(N1/2);

n2=fix(N2/2);

for i=1:N1

for j=1:N2

d=sqrt((i-n1)^2+(j-n2)^2);

h=1-exp(-d*d/(2*d0*d0));

result(i,j)=h*g(i,j);

end

result=ifftshift(result);

X2=ifft2(result);

%Fourier decentralization and inverse

final=uint8(real(X2));

%Normalize the results

figure

imshow(final);

title(['Result of GHPF with D0 = ',num2str(a)]);

end

Result and analysis:

We can find in the output figure of GLPF, the clarity increases as the value of D0 increases. Because the D0 represent the cutoff frequency. Even the noise part of the figure become much more clearly.

Then we take a closer look at the result picture and the input picture:

The left one is input figure and the right one is output figure with D0=160. By observing and comparing the details of two figures. We can find that the figure after the Gaussian filter becomes smoother. But at the edge of the figure is becoming fuzzy.

We can only know the figure become more and more darker with the value of D0 increasing. But if we have a closer look of these pictures. We will find some differences that we haven’t found.

Figure 1 Figure 2

Figure1 is the detail of output figure with D0=30

Figure2 is the detail of output figure with D0=60

Figure3 is the detail of output figure with D0=160

Figure 3

We can find that the edge of the letter ‘a’ become more and more accurate and the white band become more and more narrow.

Case3: Butterworth notch filters

Principle of Sobel filter:

Notch filters are also used to remove periodic noise. Although band stop filters can remove periodic noise, band stop filters also attenuate components other than noise. A notch filter attenuates a point, but does not lose the rest of the components.

In the spatial domain, the image has periodic noise. In the Fourier spectrum, you can see where the noise is (I've circled them in red here, they're not part of the data). If we use a band-stop filter, it is very troublesome and it will cause image loss. Notch filter can attenuate the noise directly and has good effect of denoising.

The expression is shown below. The notch filter can be obtained by shifting the center of the high-pass filter.

Here, because of the periodicity of the Fourier transform, there is no way in the Fourier spectrum that there can be a noise at a single point, there must be a noise pair symmetric about the far point. Here and The method to get the noise point is as follows:

Matlab codes:

close all;

clear all;

clc;

img = imread('Q3_3.tif');

Butterworth_11711118(img);

function Butterworth_11711118(imputImage);

f=imputImage;

f = mat2gray(f,[0 255]);

[M,N] = size(f);

P = 2*M;

Q = 2*N;

fc = zeros(M,N);

for x = 1:1:M

for y = 1:1:N

fc(x,y) = f(x,y) * (-1)^(x+y);

end

F = fft2(fc,P,Q);

H_NF = ones(P,Q);

for x = (-P/2):1:(P/2)-1

for y = (-Q/2):1:(Q/2)-1

D = 30;

v_k = 59; u_k = 77;

D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5);