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
