Activity 19 - Restoration of Blurred image

This activity will show how to restore blurred and noisy images from an a priori knowledge of degradation function and noise. The main idea of image degradation and restoration is presented using the model below:

General model in spatial domain:

where g(x,y) is the degraded image, f(x,y) is the original image, h(x,y) is the degradation function and n(x,y) is the noise model.

General model in frequency domain:

The variables presented are the Fourier transforms of the above equation.

From the above equations, we are to find the degrading function and the noise model. The noise model is described in the previous activity, so we are left with the degrading function, or the motion blurring effect. It is given by:

The image restoration procedure is summarized in the following way:

One filter that is shown in this activity is Weiner filter which is also known as the minimum mean-square error (MMSE) filter by modeling the error through statistical characterization of the noise. It minimizes the average error giving the Weiner filter:

or

where , , and K is a specified constant.

Results

The result of filtering for different types of images are presented. The first image is the original, the second is the blurred and noisy image, the last two are the two types of Weiner filter the last being specified by K = 1.

As seen, the Weiner filter restores the blurred and noisy image but it can't recover the original image itself. The phrase can still be discerned from the restoration but not fully. The last image is not a good restoration using the constant K. Even if K value is varied, no visible change in the restored image is observed.

I give myself 10 pts for understanding and finishing the activity. :)

References
M. Soriano. A19 – Restoration of blurred image. pdf file
Digital Image Processing. Gonzales and Woods.

Scilab Code

// Degradation

chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act19\');
img = gray_imread('k.bmp');
//img = gray_imread('phrase_ei.bmp');
//img = gray_imread('people1.bmp');
FTimg = fft2(img);
FT = fftshift(abs(FTimg));

a = .05;
b = .05;
T = 1;

// noise
y = grand(size(img,1), size(img,2), 'nor', 0.1, 0.1); // gauss

FTy = fft2(y);

for i = 1:size(FT,1)
for j = 1:size(FT,2)
H(i,j) = T/(%pi*(i*a + j*b)) * sin(%pi*(i*a + j*b)) * exp(-%i*(%pi*(i*a + j*b)));
end
end

G = H.*FTimg + FTy;
Sn = FTy.*conj(FTy);
Sf = FTimg.*conj(FTimg);

a = (H.^(-1));
b = ((H.*conj(H))./((H.*conj(H)) + (Sn./Sf)));
Fh1 = a.*b.*G;

K = 10000000000;
b = ((H.*conj(H))./((H.*conj(H)) + K));
Fh2 = a.*b.*G;

//scf(), imshow(FTimg, [])
scf(), imshow(abs(ifft(G)), [])
scf(), imshow(abs(ifft(Fh1)), [])
scf(), imshow(abs(ifft(Fh2)), [])
//scf(), imshow(img, [])

Activity 18 - Noise Models and basic Image Restoration

This activity allows one to digitally restore an image or by improving its appearance through a procedure of restoration and a priori knowledge or known mathematical model of degradation applied to the image. In this case, different types of noise are added on the original image.

There are different types of noise presented in this activity. These noise are experimentally obtained from the variations of temperature in the experimental area, light source flickering and also during the transmission of data through wireless connections where atmospheric disturbances may happen. Random variables comprise noise and are characterized by probability density function (PDF as discussed in Activity 4). Most common noise include:

• Gaussian
• Rayleigh
• Gamma
• Exponential
• Uniform
• Impulse (Salt and Pepper)

To restore images with these added noise, a very basic method to use is to filter in spatial domain. Filtering methods employed in this activity are as follows:

• Arithmetic
• Geometric
• Harmonic
• ContraHarmonic

where S(x,y) represent the set of coordinates in a rectangular subimage window of size m x n with center at x,y. Thus, the window S(x,y) moves through the image and proper restoration process occurs. For the ContraHarmonic Filter, Q is the order of filter with Q>0 eliminating pepper noise and Q<0 eliminating salt noise;however, it can't do both. Below is the test pattern:

RESULTS
Noise is being added to the original image. For a kind of noise added, all forms of filtering are imposed.

It can be observed that filtering the noisy images is a difficult task to do. The restored images does not contain sharp boundaries between the grayscale values. Also as seen from the PDF's of the restored images, there is really a distortion that has happened. Among the filters, the harmonic filter provides the best result except for the pepper noise. It best eliminates salt noise. The Q value used for the contraharmonic filter is 1. It better eliminates pepper than salt noise. As observed, geometric filter also does not work best for images with pepper noise.

From left to right:

original image+noise, arithmetic filtering, geometric, harmonic, contraharmonic

Gaussian noise

Rayleigh Noise



Gamma Noise



Exponential Noise



Uniform Noise



Salt Noise

Pepper noise

OTHER images

I also applied the filtering to other types of images.

Gaussian

Rayleigh

Gamma

Exponential

Uniform

Salt

Pepper

I give myself 10 pts for finishing and understanding the activity. :)

References
M. Soriano. A18 – Noise models and basic image restoration.

Scilab Code

chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act18')

function [mat] = res(arr);
mat = matrix(arr,1,-1);
endfunction

function [hist] = pdf(arr);
mx = max(arr);
hist = [];
for i = 1:mx
f = find(arr == i);
hist(i) = length(f);
end
endfunction


// constant pattern of three levels of grayscale increment
//sz = 127;
//x = -127:127;
//[X,Y] = meshgrid(x);
//s1 = (X-32)*Y;
//s2 = (X+32)*Y;
//szx = size(X,1); szy = size(Y,1); sz = szx*szy;
//s = zeros(szx, szy) + 1;
//w = 100; s(127 - w:127+w,127-w:127+w) = 128;
//w = 50; s(127 - w:127+w,127-w:127+w) = 256;
//scf(), imshow(s, [])

//
//s = gray_imread('pattern.bmp')*255 + 1;
s = gray_imread('people1.bmp')*255 + 1;
//s = gray_imread('wiki1.bmp')*255 + 1;

szx = size(s,1); szy = size(s,1); sz = szx*szy;

// noise

//Y = grand(szx, szy, 'nor', 0.1, 0.1)*255 + 1; // gauss
//Y = grand(szx, szy, 'gam', 20, 100)*255 + 1; // gamma
//Y = grand(szx, szy, 'exp', 0.1)*255 + 1; // exponential
//Y = matrix(genrayl(0.05,sz), szx, szy)*255 + 1; // rayleigh
//Y = grand(szx, szy, 'unf', 0, 0.1)*255 + 1;

//
if min(Y) < 0
Y = abs(min(Y)) + Y + 1;
end

S = s + Y;
S = round((S/max(S)) * 256);


// salt and pepper
d = 0.009;
//S = round(imnoise(s, 'salt & pepper', d));

// salt only
solt = s;
sz = size(s,1)*size(s,2);
bits = round(d*sz);
solt(round(rand(bits,1)*(sz-1)+1)) = 256;
//S = round(solt);

// pepper only
peper = s;
peper(round(rand(bits,1)*(sz-1)+1)) = 1;
S = round(peper);


// histogram
//scf(), plot(pdf(s))
//scf(), plot(pdf(S))
scf(), imshow(S, [])


// arithmetic mean filtering
filtered_ari = zeros(szx,szy);
filtered_geo = zeros(szx,szy);
filtered_har = zeros(szx,szy);
filtered_con = zeros(szx,szy);
Q = 1;
win = 3;

for i = 1:szx
for j = 1:szy
iwm = i - win;
iwp = i + win;
jwm = j - win;
jwp = j + win;
sxy = 0;
su = 0;
pr = 0;
hr = 0;
cr = 0;
fac = 1E-9;
//[i j]
if iwm <= 0 & jwm <= 0 // upper left corner
sxy = cat(2, res(S(1:i, 1:j)), res(S(i:iwp, 1:j)), res(S(1:i, j:jwp)), res(S(i:iwp,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);
//su = mean( mean(S(1:i, 1:j)) + mean(S(i:iwp, 1:j)) + mean(S(1:i, j:jwp)) + mean(S(i:iwp,j:jwp)) );

elseif iwm > 0 & iwp <= szx & jwm <= 0 // left column
sxy = cat(2, res(S(iwm:i,1:j)) , res(S(i:iwp,1:j)) , res(S(iwm:i,j:jwp)) , res(S(i:iwp,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif iwm > 0 & iwp <= szx & jwp >= szx // right column
jj = abs(j-szx);
sxy = cat(2, res(S(iwm:i,jwm:j)) , res(S(i:iwp,j:j+jj)) , res(S(iwm:i,jwm:j)) , res(S(i:iwp,j:j+jj))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif iwp > szy & jwm <= 0 // lower left corner
ii = abs(i-szy);
sxy = cat(2, res(S(iwm:i,1:j)), res(S(i:i+ii,1:j)) , res(S(iwm:i,j:jwp)) , res(S(i:i+ii,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif jwm > 0 & jwp <= szy & iwm <=0 // upper row
sxy = cat(2, res(S(1:i,jwm:j)) , res(S(i:iwp,jwm:j)) , res(S(1:i,j:jwp)) , res(S(i:iwp,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif jwm > 0 & jwp <= szy & iwp >= szy // bottom row
ii = abs(i-szy);
sxy = cat(2, res(S(iwm:i,jwm:j)) , res(S(i:i+ii,jwm:j)) , res(S(iwm:i,j:jwp)) , res(S(i:i+ii,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif iwm <= 0 & jwp > szx // upper right corner
jj = abs(j-szx);
sxy = cat(2, res(S(1:i,jwm:j)) , res(S(i:iwp,j:j+jj)) , res(S(1:i,jwm:j)) , res(S(i:iwp,j:j+jj))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);


elseif iwp > szy & jwp > szx // lower right corner
ii = abs(i-szy);
jj = abs(j-szx);
sxy = cat(2, res(S(iwm:i,jwm:j)) , res(S(i:i+ii,jwm:j)) , res(S(iwm:i,jwm:j+jj)) , res(S(i:i+ii,j:j+jj))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif i > win & i <> win & j < szx-win // no problem region
sxy = cat(2, res(S(iwm:i,jwm:j)) , res(S(i:iwp,jwm:j)) , res(S(iwm:i,j:jwp)) , res(S(i:iwp,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

end
filtered_ari(i,j) = su;
filtered_geo(i,j) = pr;
filtered_har(i,j) = hr;
filtered_con(i,j) = cr;
end
end

//scf(), imshow(s, [])
//scf(), imshow(S, [])
scf(), imshow(filtered_ari, [])
scf(), imshow(real(filtered_geo), [])
scf(), imshow(real(filtered_har), [])
scf(), imshow(real(filtered_con), [])

//scf(), plot(pdf(round(filtered_ari)))
//scf(), plot(pdf(round(filtered_geo)))
//scf(), plot(pdf(round(filtered_har)))
//scf(), plot(pdf(round(filtered_con)))
//