Fundamentals of Digital Image and Video Processing - Week 8 Solutions

Hi Coursera people,

I never knew this week would be so much interesting! I could understand almost every bit of information spoken in the class! :)

Here is the question number 7 of week 8.

In this problem, you will write a MATLAB program to compute the entropy of a given gray-scale image. Follow the instructions below to finish this problem. (1) Download the input image from here. The input is a gray-scale image with pixel values in the range [0,255]. Treat the pixel intensities in this image as symbols emitted from a DMS. (2) Build a probability model (i.e., an alphabet with associated probabilities) corresponding to this input image. Specifically, this alphabet consists of symbols {0,1,2,⋯,255}. In order to find the probabilities associated with each symbol, you will need to scan over all the pixels in this image, and for each pixel, adjust the probability associated with that pixel's intensity value accordingly, or in other words find the histogram of the image. Make sure you normalize the probability model correctly such that each probability is a real-valued number in [0,1]. (3) Compute the entropy using the formula that you have learned in class. Enter the result below to at least 2 decimal points.

Here is my attempt at the code.

A = imread('C:\~Coursera courses\Image and Video processing\Week 8\Cameraman256.bmp');
for i = 1:256
DMS(i,1) = i-1;
DMS(i,2) = 0;
end
for i = 1:256
for j = 1:256
for k = 1:256
if A(i,j) == DMS(k,1)
DMS(k,2) = DMS(k,2)+1;
end
end
end
end
sum = 0;
for i = 1:256
sum = sum + DMS(i,2);
end
for i = 1:256
prob(i) = DMS(i,2)/sum;
end
ans=0;
for i = 1:256
entropy(i) = -1 * prob(i) * log2(prob(i));
ans = ans + entropy(i);
end
ans

-Cheers,
Vijay.

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.

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.

Fundamentals of Digital Image and Video Processing - Week 5 Solutions

Hi Coursera people,

This week's videos were large and I expected the assignment also to be challenging. On the contrary, the assignment was a very simple and was a replica of week 1's assignment.

Here is the question number 7 of week 5 :

In this problem you will perform median filtering to enhance the quality of a noise corrupted image. Recall from the video lecture that median filtering is effective for removing "salt-and-pepper" noise from images. Follow the instructions below to complete this problem. (1) Download the noisy image from here. Load the noisy image into a MATLAB array and convert the type of the array from 8-bit integer 'uint8' to real number 'double'. Refer to MATLAB problems in previous homework if you need help with loading and converting images. Visualize the noisy image using the built-in MATLAB function "imshow". The function "imshow" takes as its argument either [0-255] for an 8-bit integer array (i.e., of type 'uint8'), or [0-1] for a normalized real-valued array (i.e., of type 'double'). To provide "imshow" with the correct argument, you would need either to "cast" your real-valued array into 'uint8', or normalize it by 255. (2) Perform 3x3 median filtering using the built-in MATLAB function "medfilt2". For this problem, the only argument you need to provide "medfilt2" with is the array you have created in step (1). Visualize the filtered image using "imshow". Remember to either cast the result to 'uint8' or normalize it before feeding it to "imshow". (3) Perform a second-pass median filtering on the filtered image that you have obtained from step (2). Visualize the two-pass filtered image. Compare it with the noisy input image and the 1-pass filtered image. (4) Download the noise-free image from here. Compute the PSNR values between (a) the noise-free image and the noisy input image, (b) the noise-free image and the 1-pass filtering output, and (c) the noise-free image and the 2-pass filtering output. Enter the three PSNR values in the box below. Enter the numbers to two decimal points.

And here's my attempt at the code :

A = imread('C:\xxxxxx\Week 5\digital-images-week5_quizzes-noisy.jpg');

B = im2double(A);

C = medfilt2(B,[3 3]);

D = medfilt2(C,[3 3]);

E = imread('C:\xxxxxx\Week 5\digital-images-week5_quizzes-original.jpg');

F = im2double(E);

MaxI=1;

MSE1= mean(mean((F- B).^2));

PSNR1 = 10*log10((MaxI^2)/MSE1);

MSE2 = mean(mean((F-C).^2));

PSNR2 = 10*log10((MaxI^2)/MSE2);

MSE3 = mean(mean((F-D).^2));

PSNR3 = 10*log10((MaxI^2)/MSE3);

-Cheers,

Vijay.

Fundamentals of Digital Image and Video Processing - Week 4 Solutions

Hi coursera people,

This week was very informative and very stretchy! Here is the solution to question number 8 of week 4.

Firstly, here is the question :

In this problem you will perform block matching motion estimation between two consecutive video frames. Follow the instructions below to complete this problem. (1) Download the two video frames from frame_1 and frame_2. The frames/images are of height 288 and width 352. (2) Load the frame with file name "frame_1.jpg" into a 288×352 MATLAB array using function "imread", and then convert the array type from 8-bit integer to real number using function "double" or "cast" (note that the range of intensity values after conversion is between 0 and 255). Denote by I1 the converted MATLAB array. Repeat this step for the frame with file name "frame_2.jpg" and denote the resulting MATLAB array by I2. In this problem, I2 corresponds to the current frame, and I1 corresponds to the previous frame (i.e., the reference frame). (3) Consider the 32×32 target block in I2 that has its upper-left corner at (65,81) and lower-right corner at (96,112). Note this is MATLAB coordinate convention, i.e., the first number between the parenthesis is the row index extending from 1 to 288 and the second number is the column index extending from 1 to 352. This target block is therefore a 32×32 sub-array of I2. (4) Denote the target block by Btarget. Motion estimation via block matching searches for the 32×32 sub-array of I1 that is "most similar" to Btarget. Recall in the video lectures we have introduced various forms of matching criteria, e.g., correlation coefficient, mean-squared-error (MSE), mean-absolute-error (MAE), etc. In this problem, we use MAE as the matching criterion. Given two blocks B1 and B2 both of size M×N, the MAE is defined as MAE(B1,B2)=1M×N∑Mi=1∑Nj=1|B1(i,j)−B2(i,j)|. To find the block in I1 that is most similar to Btarget in the MAE sense, you will need to scan through all the 32×32 blocks in I1, compute the MAE between each of these blocks and Btarget, and find the one that yields the smallest value of MAE. Note in practice motion search is only performed over a certain region of the reference frame, but for the sake of simplicity, we perform motion search over the entire reference frame I1 in this problem. When you find the matched block in I1, enter the following information: (1) the coordinate of the upper-left corner of the matched block in MATLAB convention. This requires two integer numbers; (2) the corresponding MAE value, which is a floating-point number. Enter the last number to two decimal points. As an example for format of answer, suppose the matched block has upper-left corner located at (1,1), and the corresponding MAE is 10.12, then you should enter 1 1 10.12 (the three numbers are separated by spaces).

And, here is my attempt at the solution :

frame_1 = imread('D:\Image processing\Week 4\digital-images-week4_quizzes-frame_1.jpg');

frame_2 = imread('D:\Image processing\Week 4\digital-images-week4_quizzes-frame_2.jpg');

I1 = double(frame_1);

I2 = double(frame_2);

Btarget = I2(65:96,81:112);

for i=1:288

if (i+31 <= 288)

for j=1:352

if (j+31 <= 352)

Btemp = I1(i:i+31,j:j+31);

err = Btarget - Btemp;

absoluteerr = abs(err);

ComputedMAE = mean2(absoluteerr);

MAEArray(i,j) = ComputedMAE;

end

end

end

end

A = min(MAEArray(:))

X = MAEArray;

[p,q] = find(X==min(X(:)))

-Cheers,

Vijay.