Fundamentals of Digital Image and Video Processing - Week 6 Solutions

Hi Coursera people,

This week's lectures were very exasperating with regards to their length! Yes, I understand that there is a time limit of 12 weeks to complete the material; but it is equally important for followers to understand what is happening right? Anyways, with what I could manage to understand, I've cooked up this week's solution. Its very surprising to learn that a few lines of code will solve the purpose...

Here is question number 6 of week 6.

This problems pertains to inverse filtering. You should review the corresponding slides in the video lectures to refresh your memory before attempting this problem. To help you understand how inverse filter is implemented and applied, we have provided you with a MATLAB script here. Download the script and the original image, and open the script using MATLAB. Once you open the script, you will see on Line 8 the statement "T = 1e-1". This defines the threshold value used in the inverse filter. The script simulates the blur due to motion and applies inverse filtering for its removal. We encourage you to try different values of the threshold and see how it affects the performance of the inverse filter. We ask you to enter the ISNR value below when the threshold is set to 0.5. Make sure you enter the number with at least 2 decimal points.

Here's the code that was given as a part of the question. I've added the code that gives the solution at the end.

% inverse filter with thresholding

clear all
close all
clc

% specify the threshold T
T = 0.5;

%% read in the original, sharp and noise-free image
original = im2double(rgb2gray((imread('C:\Image and Video processing\Week 6\original_cameraman.jpg'))));
[H, W] = size(original);

%% generate the blurred and noise-corrupted image for experiment
motion_kernel = ones(1, 9) / 9;  % 1-D motion blur
motion_freq = fft2(motion_kernel, 1024, 1024);  % frequency response of motion blur
original_freq = fft2(original, 1024, 1024);
blurred_freq = original_freq .* motion_freq;  % spectrum of blurred image
blurred = ifft2(blurred_freq);
blurred = blurred(1 : H, 1 : W);
blurred(blurred < 0) = 0;
blurred(blurred > 1) = 1;
noisy = imnoise(blurred, 'gaussian', 0, 1e-4);


%% Restoration from blurred and noise-corrupted image
% generate restoration filter in the frequency domain
inverse_freq = zeros(size(motion_freq));
inverse_freq(abs(motion_freq) < T) = 0;
inverse_freq(abs(motion_freq) >= T) = 1 ./ motion_freq(abs(motion_freq) >= T);
% spectrum of blurred and noisy-corrupted image (the input to restoration)
noisy_freq = fft2(noisy, 1024, 1024);
% restoration
restored_freq = noisy_freq .* inverse_freq;
restored = ifft2(restored_freq);
restored = restored(1 : H, 1 : W);
restored(restored < 0) = 0;
restored(restored > 1) = 1;

%% analysis of result
noisy_psnr = 10 * log10(1 / (norm(original - noisy, 'fro') ^ 2 / H / W));
restored_psnr = 10 * log10(1 / (norm(original - restored, 'fro') ^ 2 / H / W));


%% visualization
figure; imshow(original, 'border', 'tight');
figure; imshow(blurred, 'border', 'tight');
figure; imshow(noisy, 'border', 'tight');
figure; imshow(restored, 'border', 'tight');
figure; plot(abs(fftshift(motion_freq(1, :)))); title('spectrum of motion blur'); xlim([0 1024]);
figure; plot(abs(fftshift(inverse_freq(1, :)))); title('spectrum of inverse filter'); xlim([0 1024]);

%% Calculation of ISNR

e1=original-noisy;
e2=original-restored;
E1=mean2(e1.*e1);
E2=mean2(e2.*e2);
result=10*log(E1/E2)/log(10)




P.S. : The answer is varying with regards to the second decimal as this code is run several times on the same machine. Don't ask me why! I'm as clueless as you are :p
Just type "2.85" in the answer area, and you get a full 3 points!

-Cheers,
Vijay.