Project 1: Image Filtering and Hybrid Images

This script show cases self implemented filtering functions, and uses them to create example cases of hybrid images.

Part 1 - Image Filtering

Here I create 2 filtering functions that behaves exactly the same way as the MATLAB built in function imfilter. One implementation is in real space, the other in fourier space.

Real Space Filtering

Filtering (or Convolution) can be written as an inner product of a partial Circulant matrix and a filter vector . This code exploits this formulation by calculating individual rows of (which correspond to shifts of the image) and multiplying them with the scalar element of the filter that corresponds to that shift.

function output = my_imfilter(image, filter)
% Real space image filtering.
% Only accept odd filters.
if sum(mod(size(filter),2)) ~= 2 
    error('Filter must have odd dimensions.');
end 
% Make sure inputs are double.
image = double(image);
filter = double(filter);

% Transpose filter to accomodate MATLAB x,y vs i,j notation.
filter = filter';

% Get input sizes explicitly.
imagexsize = size(image,1);
imageysize = size(image,2);
filterxspan = (size(filter,1)-1)/2;
filteryspan = (size(filter,2)-1)/2;

% Boundary conditions - Zero padded boundary.
image = padarray(image,[filterxspan filteryspan],0,'both');

% Pre-allocate output.
output = zeros(size(image));

% Get shift indices.
[xgrid, ygrid] = meshgrid(-filterxspan:filterxspan,-filteryspan:filteryspan);

% Loop over filter elements.
for i = 1:numel(filter)
    
% Shift image based on the index of the current filter element.
    current = circshift(image,[-xgrid(i) -ygrid(i)]);
    
% Multiply with current filter element.
    current = current.*filter(i);
    
% Sum over all filter elements to get the filtered image.
    output = output + current;
    
end
% Resize output to original input dimensions.
output = output(filterxspan+1:filterxspan+imagexsize,filteryspan+1:filteryspan+imageysize,:);

Fourier Space Filtering

Any circulant matrix can be diagonalized by the Fourier matrix F as:

This fact can be exploited to achieve multiplication of a vector with in O(nlogn) time using FFT. The following scheme can be used for this purpose, where is the first column and is a vector containing the main diagonal of :

where is element-wise multiplication. Notice the above algortihm is actually the Convolution Theorem in disguise, where is the Fourier Transform operator:

The mirror mismatch between convolution and correlation operations can be accounted for by conjugation of the filter in Fourier space.

function output = my_imfilter_fft(image, filter)
% Image filtering using Convolution Theorem.
% Only accept odd filters.
if sum(mod(size(filter),2)) ~= 2 
    error('Filter must have odd dimensions.');
end
% Make sure inputs are double.
image = double(image);
filter = double(filter);

% Get input sizes explicitly.
imagexsize = size(image,1);
imageysize = size(image,2);
filterxspan = (size(filter,1)-1)/2;
filteryspan = (size(filter,2)-1)/2;

% Boundary conditions - Zero padded boundary.
image = padarray(image,[filterxspan filteryspan],0,'both');

% Pad filter to image dimensions for the FFT method.
filter = padarray(filter,[size(image,1) size(image,2) size(image,3)]-[size(filter) size(filter,3)],0,'post'); 

% Shift filter to FFT indexing.
filter = circshift(filter,[-filterxspan -filteryspan]);

% Convolution Theorem
output = ifftn(fftn(image).*conj(fftn(filter)));

% Resize output to original input dimensions.
output = output(filterxspan+1:filterxspan+imagexsize,filteryspan+1:filteryspan+imageysize,:);

Testing Accuracy and Performance

Let us load a test image first.

test_image = im2single(imread('../data/cat.bmp'));
test_image = imresize(test_image, 0.7, 'bilinear');% Resize for faster run.

figure
imshow(test_image);

First lets us see if our functions are accurate by comparing them with the output from imfilter.

random_filter = rand(5);
output1 = imfilter(test_image, random_filter);
output2 = my_imfilter(test_image, random_filter);
output3 = my_imfilter_fft(test_image, random_filter);

realspaceerror = norm(output2(:)-output1(:))/norm(output1(:));
fourierspaceerror = norm(output3(:)-output1(:))/norm(output1(:));

fprintf('Real Space Error = %f \n',realspaceerror);
fprintf('Fourier Space Error = %f \n',fourierspaceerror);
Real Space Error = 0.000000 
Fourier Space Error = 0.000000 

Now we can compare the speed of the 3 functions on a big image. The church image is 20 megapixels. Keep in mind that while the performance of MATLAB built-in function and real space implementation greatly deteriorates with larger filters, the Fourier space implementation only depends on the image size.

large_image = imread('../data/church.jpg');
random_filter = rand(5);

tic;
output1 = imfilter(large_image, random_filter);
t1 = toc;

tic;
output2 = my_imfilter(large_image, random_filter);
t2 = toc;

tic;
output3 = my_imfilter_fft(large_image, random_filter);
t3 = toc;

figure
imshow(large_image);

figure
bar([t1,t2,t3]);
set(gca,'xticklabel',{'MATLAB Built-in','Real Space','Fourier Space'});
title('Comparison of Runtimes(s)');
Warning: Image is too big to fit on screen; displaying at 17% 

Example Filtering Cases

Identity Filter

identity_filter = [0 0 0; 0 1 0; 0 0 0];
identity_image  = my_imfilter(test_image, identity_filter);

figure
imshow(identity_image);

Small Blur Filter

blur_filter = [1 1 1; 1 1 1; 1 1 1];
blur_filter = blur_filter / sum(sum(blur_filter));
blur_image = my_imfilter(test_image, blur_filter);

figure
imshow(blur_image);

Large Blur Filter - Exploit Seperability

large_1d_blur_filter = fspecial('Gaussian', [25 1], 10);
large_blur_image = my_imfilter(test_image, large_1d_blur_filter);
large_blur_image = my_imfilter(large_blur_image, large_1d_blur_filter');

figure
imshow(large_blur_image);

Sobel Filter

sobel_filter = [-1 0 1; -2 0 2; -1 0 1];
sobel_image = my_imfilter(test_image, sobel_filter);

figure
imshow(sobel_image + 0.5);

High Pass Filter

laplacian_filter = [0 1 0; 1 -4 1; 0 1 0];
laplacian_image = my_imfilter(test_image, laplacian_filter);

figure
imshow(laplacian_image + 0.5);

Part 2 - Hybrid Images

Hybrid images are created by melding the high frequency component of one image with the low frequency component of another. The images need to be aligned well for the product to be attractive.

image1 = im2single(imread('../data/dog.bmp'));
image2 = im2single(imread('../data/cat.bmp'));

Blur filter for low frequencies.

cutoff_frequency = 7;
filter = fspecial('Gaussian', cutoff_frequency*4+1, cutoff_frequency);

Low Frequency Component of the First Image

low_frequencies = my_imfilter(image1,filter);

figure
imshow(low_frequencies);

High Frequency Component of the Second Image

low_frequencies2 = my_imfilter(image2,filter);
high_frequencies = image2 - low_frequencies2;

figure
imshow(high_frequencies + 0.5);

Hybrid Image

hybrid_image = high_frequencies + low_frequencies;
vis = vis_hybrid_image(hybrid_image);
figure
imshow(vis);


Published with MATLABĀ® R2016b