Contents

Implement your own code for image filtering, my_imfilter() in MATLAB or C++ (You will get extra credit if you use C++). It must (1) support grayscale and color images (2) support arbitrary shaped filters, as long as both dimensions are odd (e.g. 7x9 filters but not 4x5 filters) (3) pad the input image with zeros or reflected image content and (4) return a filtered image which is the same resolution as the input image. Test your program on the provided images.

Forbidden MATAB functions which you cannot use in your final code: imfilter(), filter2(), conv2(), nlfilter(), colfilt().

ccc;

1.1. Image filtering

simple image filtering using 'my_imfilter'

img_idx = 1;
img_full_path = get_img_path(img_idx);
fprintf('\nImage filtering [%s] LOADED. \n', img_full_path);
raw_img = imread(img_full_path);
fw = 9;
fh = 9;
filter = ones(fh, fw);
filter = filter/sum(sum(filter));
fprintf('Evaluating filtering...'); tic;
fimg_matlab = imfilter(raw_img, filter); % <== this is just for comparison
fimg_mine = my_imfilter(raw_img, filter);
fprintf('finished. %.e sec \n', toc);

fig = figure(1); clf;
subaxes(fig, 2, 2, 1);
imshow(raw_img);
title('Raw Image', 'FontSize', 11);
subaxes(fig, 2, 2, 2);
imshow(fimg_mine);
title('Filtered Image (mine)', 'FontSize', 11);
subaxes(fig, 2, 2, 3);
imshow(fimg_matlab);
title('Filtered Image (solution from MATLAB)', 'FontSize', 11);
subaxes(fig, 2, 2, 4);
diff_img = fimg_matlab - fimg_mine;
imshow(diff_img);
% compute error
nr_ch = size(diff_img, 3);
diff = 0;
for i = 1:nr_ch
    temp_img = diff_img(:, :, i);
    diff = diff + norm(double(temp_img));
end
title(['Error btw Mine and Matlab is ' num2str(diff, '%.3e')], 'FontSize', 11);
drawnow;

Image filtering [../../images/bike/bicycle.bmp] LOADED. 
Evaluating filtering...finished. 2e+00 sec 

1.2. Hybrid image

Generate hybrid image

img_idx1 = 1;
img_idx2 = 2;
img_full_path1 = get_img_path(img_idx1);
img_full_path2 = get_img_path(img_idx2);

img1 = imread(img_full_path1);
img2 = imread(img_full_path2);
fprintf('\nHybrid image [%s] + [%s] ARE LOADED. \n', img_full_path1, img_full_path2);

% 2. Do hybrid filtering using 'img1' and 'img2'
%
% We will use [0~1] converted double image
%
img1_double = im2double(img1);
img2_double = im2double(img2);

% Image for low pass and high pass
img4lp_double = img2_double;
img4hp_double = img1_double;

lp_len = 21; lp_dev = 5;
lp_filter = (fspecial('gaussian', [lp_len, lp_len], lp_dev));
hp_len = 9; hp_dev = 1;
hp_half = (hp_len-1)/2;
hp_filter = padarray(1, [hp_half hp_half]) - (fspecial('gaussian', [hp_len, hp_len], hp_dev));

sum_lp = sum(sum(lp_filter));
sum_hp = sum(sum(hp_filter));
% fprintf('sum_lp: %.2e / sum_hp: %.2e \n', sum_lp, sum_hp);

plot_filter = 0;
if plot_filter
    fig = figure(2); clf;
    subaxes(fig, 1,2,1);
    imagesc(lp_filter); colorbar; axis square;
    title('Low pass filter');
    subaxes(fig, 1,2,2);
    imagesc(hp_filter); colorbar; axis square;
    title('High pass filter');
end

tic;
fprintf('Evaluating filtering...'); tic;
lp_img_double = my_imfilter(img4lp_double, lp_filter);
hp_img_double = my_imfilter(img4hp_double, hp_filter);
fprintf('finished. %.e sec \n', toc);
hyb_img_double = lp_img_double + hp_img_double;
resize_rate = 0.75;
resize_rate1 = resize_rate;
r_hyb_img_double1 = imresize(hyb_img_double, resize_rate1, 'bicubic');
resize_rate2 = resize_rate^2;
r_hyb_img_double2 = imresize(hyb_img_double, resize_rate2, 'bicubic');
resize_rate3 = resize_rate^3;
r_hyb_img_double3 = imresize(hyb_img_double, resize_rate3, 'bicubic');
resize_rate4 = resize_rate^4;
r_hyb_img_double4 = imresize(hyb_img_double, resize_rate4, 'bicubic');
resize_rate5 = resize_rate^5;
r_hyb_img_double5 = imresize(hyb_img_double, resize_rate5, 'bicubic');

% plot
fig = figure(2); clf; scrsz = get(0,'ScreenSize');
xstart = 0.1; ystart = 0.3;
xsize  = 0.5; ysize  = 0.5;
set(fig,'Position', [xstart*scrsz(3), ystart*scrsz(4), xsize*scrsz(3), ysize*scrsz(4)]);
subaxes(fig, 2,3,1);
imshow(img4lp_double);
title('Original image for Low-pass', 'FontSize', 15);
subaxes(fig, 2,3,2);
imshow(img4hp_double);
title('Original image for High-pass', 'FontSize', 15);
subaxes(fig, 2,3,3);
imshow(lp_img_double);
title('Low passed image', 'FontSize', 15);
subaxes(fig, 2,3,4);
imshow(hp_img_double+0.5);
title('High passed image + 0.5', 'FontSize', 15);
subaxes(fig, 2,3,5);
imshow(hyb_img_double);
title('Hybrid image', 'FontSize', 15);
drawnow;

% plot image pyramid
clear imagesCellArray
mand = hyb_img_double;
dim_r = 1;
dim_c = 5+1;
[imagesCellArray{1:dim_r,1:dim_c}] = deal(mand); % create smaller images by imresize
str = cell(dim_c);
for iRow = 1:dim_r
    for iCol = 1:dim_c
        r_rate = resize_rate^(iCol-1);
        imagesCellArray{iRow,iCol} = imresize(imagesCellArray{iRow,iCol}, r_rate, 'bicubic');
        if iCol == 1
            str{iCol} = sprintf('Original image');
        else
            str{iCol} = sprintf('rate: %.2f', r_rate);
        end
    end
end

% plot with imshowTruesize - true aspect ratio is preserved
margins = [25 25];
Handles = imshowTruesize(imagesCellArray,margins);

for iRow = 1:dim_r
    for iCol = 1:dim_c
        % axis(Handles.hSubplot(iRow,iCol), 'on')
        title(Handles.hSubplot(iRow,iCol), sprintf(str{iCol}));
    end
end

Hybrid image [../../images/bike/bicycle.bmp] + [../../images/bike/motorcycle.bmp] ARE LOADED. 
Evaluating filtering...finished. 3e+00 sec 

1.3 Frequency analysis

For your favorite result, you should also illustrate the process through frequency analysis. Show the log magnitude of the Fourier transform of the two input images, the filtered images, and the hybrid image. In MATLAB, you can compute and display the 2D Fourier transform with: imagesc(log(abs(fftshift(fft2(gray_image)))))

g_img4lp_double = rgb2gray(img4lp_double);
g_img4hp_double = rgb2gray(img4hp_double);
g_hp_img_double = rgb2gray(lp_img_double);
g_lp_img_double = rgb2gray(hp_img_double);
g_hyb_img_double = rgb2gray(hyb_img_double);
g_r_hyb_img_double1 = rgb2gray(r_hyb_img_double1);
g_r_hyb_img_double2 = rgb2gray(r_hyb_img_double2);
g_r_hyb_img_double3 = rgb2gray(r_hyb_img_double3);
g_r_hyb_img_double4 = rgb2gray(r_hyb_img_double4);
g_r_hyb_img_double5 = rgb2gray(r_hyb_img_double5);

fig = figure(4); clf; scrsz = get(0,'ScreenSize');
xstart = 0.1; ystart = 0.3;
xsize  = 0.5; ysize  = 0.5;
set(fig,'Position', [xstart*scrsz(3), ystart*scrsz(4), xsize*scrsz(3), ysize*scrsz(4)]);
subaxes(fig, 3,4,1);
imagesc(log(abs(fftshift(fft2(g_img4lp_double))))); axis square;
title('Log magnitude of Original image4LP', 'FontSize', 12);
subaxes(fig, 3,4,2);
imagesc(log(abs(fftshift(fft2(g_img4hp_double))))); axis square;
title('Log magnitude of Original image4HP', 'FontSize', 12);
subaxes(fig, 3,4,3);
imagesc(log(abs(fftshift(fft2(g_hp_img_double))))); axis square;
title('Log magnitude of Low pass image', 'FontSize', 12);
subaxes(fig, 3,4,4);
imagesc(log(abs(fftshift(fft2(g_lp_img_double))))); axis square;
title('Log magnitude of High pass image', 'FontSize', 12);
subaxes(fig, 3,4,5);
imagesc(log(abs(fftshift(fft2(g_hyb_img_double))))); axis square;
title('Log magnitude of hybrid image', 'FontSize', 12);
subaxes(fig, 3,4,6);
imagesc(log(abs(fftshift(fft2(g_r_hyb_img_double1))))); axis square;
title(sprintf('Log magnitude of scaled down (%.2f)', resize_rate1), 'FontSize', 12);
subaxes(fig, 3,4,7);
imagesc(log(abs(fftshift(fft2(g_r_hyb_img_double2))))); axis square;
title(sprintf('Log magnitude of scaled down (%.2f)', resize_rate2), 'FontSize', 12);
subaxes(fig, 3,4,8);
imagesc(log(abs(fftshift(fft2(g_r_hyb_img_double3))))); axis square;
title(sprintf('Log magnitude of scaled down (%.2f)', resize_rate3), 'FontSize', 12);
subaxes(fig, 3,4,9);
imagesc(log(abs(fftshift(fft2(g_r_hyb_img_double4))))); axis square;
title(sprintf('Log magnitude of scaled down (%.2f)', resize_rate4), 'FontSize', 12);
subaxes(fig, 3,4,10);
imagesc(log(abs(fftshift(fft2(g_r_hyb_img_double5))))); axis square;
title(sprintf('Log magnitude of scaled down (%.2f)', resize_rate5), 'FontSize', 12);

drawnow;

2.1 Zero-crossing of LoG

Find edge by finding zero-crossing of LoG

img_idx = 2; % 1, 2, 3, 4, 5, 6
img_full_path = get_img_path(img_idx);
img_rgb    = imread(img_full_path);
img_double = im2double(img_rgb);
fprintf('\nZero-crossing of LoG [%s] LOADED. \n', img_full_path);

% 2. LoG Filter
fprintf('Evaluating LoG..'); tic;
% log_len = 9;
% log_dev = 1.4;
% log_filter = (fspecial('log', [log_len, log_len], log_dev));
% log_img_double = my_imfilter(img_double, log_filter);
log_img_double = my_imfilter4edge(img_double, 'log');
fprintf('Finished, took %.1f sec \n', toc);

% 3. Find edge
edge_param = get_edge_param(img_idx);
edge_img = find_edge(log_img_double, edge_param);

% plot
fig = figure(5); clf; scrsz = get(0,'ScreenSize');
xstart = 0.1; ystart = 0.3;
xsize  = 0.7; ysize  = 0.4;
set(fig,'Position', [xstart*scrsz(3), ystart*scrsz(4), xsize*scrsz(3), ysize*scrsz(4)]);
subaxes(fig, 1,3,1);
imshow(img_double);
subaxes(fig, 1,3,2);
imshow(log_img_double+0.5);
title('LoG filtered image with 0.5 bias');
subaxes(fig, 1,3,3);
imshow(edge_img);
title('Black dots indicate edge');
drawnow;

Zero-crossing of LoG [../../images/bike/motorcycle.bmp] LOADED. 
Evaluating LoG..Finished, took 1.3 sec 

2.2 Zero-crossing of LOBG

- And the zero-crossings of the LOBG (Laplacian of Bilateral Gaussian) - You can derive the LOBG equation and design the kernel similarly to LOG case ? Write a single algorithm that has two kernel options, LOG and LOGB. ? Test your algorithm on the test images by varying the parameters of each algorithm (Try 9x9 kernel size with Gaussian std = 1.4).

img_idx = 2; % 1, 2, 3, 4, 5, 6
img_full_path = get_img_path(img_idx);
img_rgb    = imread(img_full_path);
img_double = im2double(img_rgb);
fprintf('\nZero-crossing of LoBG [%s] LOADED. \n', img_full_path);

% 2. BiLoG Filter
fprintf('Evaluating LoBG..'); tic;
% gauss_len = 9;
% gauss_dev = 1.4;
% gauss_filter = (fspecial('gaussian', [gauss_len, gauss_len], gauss_dev));
% blog_img_double = my_imbifilter(img_double, gauss_filter);
blog_img_double = my_imfilter4edge(img_double, 'lobg');
fprintf('Finished, took %.1f sec \n', toc);

% 3. Find edge
edge_param = get_edge_param(img_idx);
edge_img = find_edge(blog_img_double, edge_param);

% plot
fig = figure(6); clf; scrsz = get(0,'ScreenSize');
xstart = 0.2; ystart = 0.3;
xsize  = 0.7; ysize  = 0.4;
set(fig,'Position', [xstart*scrsz(3), ystart*scrsz(4), xsize*scrsz(3), ysize*scrsz(4)]);
subaxes(fig, 1,3,1);
imshow(img_double);
subaxes(fig, 1,3,2);
imshow(blog_img_double+0.5);
title('LoBG filtered image with 0.5 bias');
subaxes(fig, 1,3,3);
imshow(edge_img);
title('Zero-crossing of LoBG Black dots indicate edge');
drawnow;

Zero-crossing of LoBG [../../images/bike/motorcycle.bmp] LOADED. 
Evaluating LoBG..Finished, took 142.8 sec 

2.3 Add error to image and compare LoG and LoBG

? Add Gaussian noise with std= 20 to the test images and repeat the test on them. ? Show and compare your results, and discuss on them. Give any ideas for improving the performance.

img_idx = 2; % 1, 2, 3, 4, 5, 6
img_full_path = get_img_path(img_idx);
img_rgb    = imread(img_full_path);
img_noise_rgb = uint8(double(img_rgb) + 20*randn(size(img_rgb)));
% img_noise_rgb = imnoise(img_rgb, 'Gaussian', 0, 20);
img_double = im2double(img_rgb);
img_noise_double = im2double(img_noise_rgb);
fprintf('\nComparison of LoG and LoBG [%s] LOADED. \n', img_full_path);


% 1. Zero crossing of LoG Filter
fprintf('Evaluating LoG..'); tic;
% log_len = 9;
% log_dev = 1.4;
% log_filter = (fspecial('log', [log_len, log_len], log_dev));
% log_img_double = my_imfilter(img_noise_double, log_filter);
log_img_double = my_imfilter4edge(img_noise_double, 'log');
fprintf('Finished, took %.1f sec \n', toc);
edge_param = get_edge_param(img_idx);
edge_log_img = find_edge(log_img_double, edge_param);


% 2. Zero crossing of LoBG Filter
fprintf('Evaluating LoBG..'); tic;
% gauss_len = 9;
% gauss_dev = 1.4;
% gauss_filter = (fspecial('gaussian', [gauss_len, gauss_len], gauss_dev));
% blog_img_double = my_imbifilter(img_noise_double, gauss_filter);
blog_img_double = my_imfilter4edge(img_noise_double, 'lobg');
fprintf('Finished, took %.1f sec \n', toc);
edge_param = get_edge_param(img_idx);
edge_lobg_img = find_edge(blog_img_double, edge_param);

% 3. Zero crossing of BiLoG Filter (with pre-processing)
pre_gauss_len = 3;
pre_gauss_dev = .5; % 1
pre_gauss_filter = (fspecial('gaussian', [pre_gauss_len, pre_gauss_len], pre_gauss_dev));
fprintf('Evaluating LoBG with preprocess..'); tic;
img_double_pre_gauss = my_imfilter(img_noise_double, pre_gauss_filter);
% gauss_len = 9;
% gauss_dev = 1.4;
% gauss_filter = (fspecial('gaussian', [gauss_len, gauss_len], gauss_dev));
% blog_pre_img_double = my_imbifilter(img_double_pre_gauss, gauss_filter);
blog_pre_img_double = my_imfilter4edge(img_double_pre_gauss, 'lobg');
fprintf('Finished, took %.1f sec \n', toc);
edge_param = get_edge_param(img_idx);
edge_lobg_pre_img = find_edge(blog_pre_img_double, edge_param);

% plot
fig = figure(7); clf; scrsz = get(0,'ScreenSize');
xstart = 0.2; ystart = 0.1;
xsize  = 0.6; ysize  = 0.7;
set(fig,'Position', [xstart*scrsz(3), ystart*scrsz(4), xsize*scrsz(3), ysize*scrsz(4)]);
subaxes(fig, 2,3,1);
imshow(img_double);
title('Original image');
subaxes(fig, 2,3,2);
imshow(img_noise_double);
title('Noisy image');
subaxes(fig, 2,3,3);
imshow(edge_log_img);
title('LoG edge');
subaxes(fig, 2,3,4);
imshow(edge_lobg_img);
title('LoBG edge');
subaxes(fig, 2,3,5);
imshow(edge_lobg_pre_img);
title('Preprocessed LoBG edge');
drawnow;

Comparison of LoG and LoBG [../../images/bike/motorcycle.bmp] LOADED. 
Evaluating LoG..Finished, took 1.4 sec 
Evaluating LoBG..Finished, took 142.7 sec 
Evaluating LoBG with preprocess..Finished, took 161.9 sec 

Published with MATLAB® R2014b