Lab 4：Image Restoration
In this lab, I review the principle of several kinds of filters. And use the histograms of different figure to determine which kind of noise has been added on the figure. Then I try to use different methods to reduce the noise to restore the original figure.Introduction: Question1: Types of noise:
By obverting the histograms of four input figures, we can determine the noise types of them.
We can know the first figure has pepper noise, the second input figure has salt noise, for first two figure, we can use median filter to reduce the noise. The third input figure has alt-and-pepper noise, we can use adaptive median filter. The fourth input figure has Gaussian noise, we can use Gaussian filter to restore the figure.
Matlab codes:
function f = RAMF(img) % adaptive median filter
[Im,In] = size(img);
nmin = 3;
nmax = 9;
Imf = img;
I_ex = [zeros((nmax-1)/2,In+(nmax-1));zeros(Im,(nmax-1)/2),img,zeros(Im,(nmax-1)/2);zeros((nmax-1)/2,In+(nmax-1))];
for x = 1:Im
for y = 1:In
for n = nmin:2:nmax
Sxy = I_ex(x+(nmax-1)/2-(n-1)/2:x+(nmax-1)/2+(n-1)/2,y+(nmax-1)/2-(n-1)/2:y+(nmax-1)/2+(n-1)/2);
Smax = max(max(Sxy));% Find the maximum number of pixels in the window
Smin = min(min(Sxy));%Find the minimum number of pixels in the window
Smed = median(median(Sxy));%Find the middle value of pixels in the window
% Determine if the median is a noise point
if Smed > Smin && Smed < Smax
% if the median is both greater than the minimum and less than the maximum
% If is, exit the if statement, increase the window size, and judge again
% is not, then the original value of the point is not a noise point
if Imf(x,y) <= Smin || Imf(x,y) >= Smax
% if the original value of the point is both greater than the minimum and less than the maximum, then it is not
%If not, then the output value, that is, not for processing
%If so, the median is printed
Imf(x,y) = Smed;
end
break
end
end
Imf(x,y) = Smed;
end
end
f = Imf;
end
function L=GaussianLow_11711118(I,setD0)
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));
end
Result and analysis:
We can see that the output figure had lost some high frequency information. And another question is the figure had an edge which I cannot reduce. For instance, like on the top edge of the output figure of the third figure has a line:
Question2: full inverse filtering Principle of full inverse filtering:
We can find that the input figure of question 2 is a figure blurred by atmosphere turbulence. So that we can use full inverse filter to restore the figure.
Matlab codes:
image_d=imread('Q4_2.tif');
ff=im2double(image_d);%The image gray value is normalized to 0-1
% Fourier transfrom
f_Id=fft2(ff);
f_Id=fftshift(f_Id);
fH_Id=f_Id;
[M,N]=size(fH_Id);
% set the threshold
threshold=90;
if threshold>M/2
fH_Id=fH_Id./(H+eps);
else
%Filter within a certain radius
for i=1:M
for j=1:N
if sqrt((i-M/2).^2+(j-N/2).^2)<threshold
fH_Id(i,j)=fH_Id(i,j)./(H(i,j)+eps);
end
end
end
end
% Perform the inverse Fourier transform
fH_Id1ftshift(fH_Id);
f_new=ifft2(fH_Id1);
f_new=uint8(abs(f_new)*255);
imshow(f_new);
title('output figure with threshold=90');
Result and analysis:
With the threshold value is 90, we can find that the effect of atmosphere turbulence had been reduced.
Question3: Butterworth notch filtersWe can see that the first two picture have noise easily. Then we can take a look at the histogram and frequency domain to determine which kind of noise the input figures have.
We can find that the frequency domain figure can give us nothing. So we need to look at the histogram, we can guess the noise is Gaussian noise.
So, we use the Gaussian lowpass filter to reduce the noise.
This figure is the first two figure pass the Gaussian filter. We can see that the noise has been reduced.
Matlab codes:
clc;
clear;
i1=imread('Q4_3_1.tiff');
i2=imread('Q4_3_2.tiff');
i3=imread('Q4_3_3.tiff');
o1=fftshift(fft2(i1));
o2=fftshift(fft2(i2));
o3=fftshift(fft2(i3));
figure(1)
subplot(231),imshow(i1);
subplot(234),imshow(o1);
subplot(232),imshow(i2);
subplot(235),imshow(o2);
subplot(233),imshow(i3);
subplot(236),imshow(o3);
figure(4)
subplot(221)
imhist(i1);
subplot(222)
imhist(i2);
subplot(223)
imhist(i3);
n1 = RAMF_11711118(i1);
n2 = RAMF_11711118(i2);
figure(5)
subplot(121)
imhist(n1);
subplot(122)
imhist(n2);
figure(6)
subplot(221),imshow(i1);
subplot(222),imshow(n1);
subplot(223),imshow(i2);
subplot(224),imshow(n2);
G1 = GaussianLow_11711118(i1,160);
G2 = GaussianLow_11711118(i2,160);
figure(7)
subplot(121)
imhist(G1);
subplot(122)
imhist(G2);
figure(8)
subplot(221),imshow(i1);
subplot(222),imshow(G1);
subplot(223),imshow(i2);
subplot(224),imshow(G2);
Answer and analysis of questions:
We can take at the question one:
Because the pepper noise is showed as small black points in the figure and slat noise is showed as white points. By the formula we get:
We need to make the black points become brighter; the white points become darker. So we need to thins the dark part for Q>0, and thickens the black part for Q<0?
We can take at the question two:
We can take a look at the formula:
If the σn2 is larger than the actual global variance, the value of fx,y will be decreasing. So the figure may become darker. And become brighter when σn2 is small than the actual global variance
