[数字图像处理]Image Restoration实验报告

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.

[数字图像处理]Image Restoration实验报告

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

Imf(x,y) = Smed;

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:

[数字图像处理]Image Restoration实验报告

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:

[数字图像处理]Image Restoration实验报告

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

% 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:

[数字图像处理]Image Restoration实验报告

With the threshold value is 90, we can find that the effect of atmosphere turbulence had been reduced.

Question3: Butterworth notch filters

We 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.

[数字图像处理]Image Restoration实验报告

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.

[数字图像处理]Image Restoration实验报告

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:

[数字图像处理]Image Restoration实验报告

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:

[数字图像处理]Image Restoration实验报告

If the σn2 [数字图像处理]Image Restoration实验报告 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