Image filtering + Hybrid image

Goal

: 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().

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;

Result

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

Result

: 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)))))

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;

Result 

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).

Result

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;

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;

Result

Full code

hw1.zip