Fundamentals of Digital Image and Video Processing - Week 7 Solutions

Hi Coursera people,

After hours of struggle, I had found the solution to question number 7 of week 7! There is absolutely no material available online regarding calculation of frequency response of the CLS filter! Finally, the answer sprung up from the discussions forum of the course itself :) I must thank my fellow Courserans for helping me out with the code. Here goes the question number 7 of week 7.

In this problem, you will implement the Constrained Least Squares (CLS) filter and examine its performance when the regularization parameter is set at different values. You will be provided with the original image and a set of MATLAB files. Follow the instructions below to finish this problem. (1) Download the original image and the MATLAB code from here. Place the original image and all the provided MATLAB files in the same directory. (2) The file "wrapper.m" is the entry or the "main" code. It loads the original image, applies a motion blur to it, and degrades the image by adding noise. The 17th line in "wrapper.m" sets the value of the regularization parameter "alpha". (3) The MATLAB file "cls_restoration.m" has an incomplete implementation of the CLS filter. You need to un-comment line 24 in "cls_restoration.m" and complete the implementation of the CLS filter. (4) After you complete the implementation of the CLS filter, you should run "wrapper.m" with different values of alpha. Specifically, we ask you to try the following values of alpha: {0.0001, 0.001, 0.01, 0.1, 1, 10, 100}. For each value of alpha, we ask you to compute the Improvement in SNR (ISNR). Note that the computation of ISNR involves there images: the original image, the blurred and noisy image, and the restored image. After you obtain the ISNR values, enter in the box below the largest ISNR value. Enter the number with at least two decimal points.

Here is my attempt at the code :

wrapper.m

clear all
close all

%% Simulate 1-D blur and noise
image_original = im2double(imread('C:\Image and Video processing\Week 7\downloaded codes\Cameraman256.bmp', 'bmp'));
[H, W] = size(image_original);
blur_impulse = fspecial('motion', 7, 0);
image_blurred = imfilter(image_original, blur_impulse, 'conv', 'circular');
noise_power = 1e-4;
randn('seed', 1);
noise = sqrt(noise_power) * randn(H, W);
image_noisy = image_blurred + noise;

figure; imshow(image_original, 'border', 'tight');
figure; imshow(image_blurred, 'border', 'tight');
figure; imshow(image_noisy, 'border', 'tight');

%% CLS restoration
alpha = 0.0001;  % you should try different values of alpha
image_cls_restored = cls_restoration(image_noisy, blur_impulse, alpha);
figure; imshow(image_cls_restored, 'border', 'tight');

%% computation of ISNR

e1=image_original-image_noisy;
e2=image_original-image_cls_restored;
E1=mean2(e1.*e1);
E2=mean2(e2.*e2);
result=10*log(E1/E2)/log(10)


cls_restoration.m

function image_restored = cls_restoration(image_noisy, psf, alpha)

%% find proper dimension for frequency-domain processing
[image_height, image_width] = size(image_noisy);
[psf_height, psf_width] = size(psf);
dim = max([image_width, image_height, psf_width, psf_height]);
dim = next2pow(dim);

%% frequency-domain representation of degradation
psf = padarray(psf, [dim - psf_height, dim - psf_width], 'post');
psf = circshift(psf, [-(psf_height - 1) / 2, -(psf_width - 1) / 2]);
H = fft2(psf, dim, dim);

%% frequency-domain representation of Laplace operator
Laplace = [0, -0.25, 0; -0.25, 1, -0.25; 0, -0.25, 0];
Laplace = padarray(Laplace, [dim - 3, dim - 3], 'post');
Laplace = circshift(Laplace, [-1, -1]);
C = fft2(Laplace, dim, dim);

%% Frequency response of the CLS filter
% Refer to the lecture for frequency response of CLS filter
% Complete the implementation of the CLS filter by uncommenting the
% following line and adding appropriate content

R = conj(H)./(abs((H.*H))+(alpha*abs((C.*C))));

%% CLS filtering
Y = fft2(image_noisy, dim, dim);
image_restored_frequency = R .* Y;
image_restored = ifft2(image_restored_frequency);
image_restored = image_restored(1 : image_height, 1 : image_width);


next2pow.m

function result = next2pow(input)
if input <= 0
    fprintf('Error: input must be positive!\n');
    result = -1;
else
    index = 0;
    while 2 ^ index < input
        index = index + 1;
    end
    result = 2 ^ index;
end


And don't get freaked if you get negative ISNR values as the result; it is absolutely normal. Here are the observations for the various values of alpha.
0.0001    -5.9191
0.001      -1.5292
0.01        3.4933
0.1          4.3048
1             2.1471
10           0.4866

-Cheers,
Vijay.