Sunday, December 20, 2015

Fourier in image processing (I)

before indulge into the frequency field of image processing, two points must be manifested:

1. the mathematical rigor is a good thing, however in practice just stick to it in some extent in case of seeking trouble for oneself. so don't be panic if you see someone exchange summation and integration freely.

2. casualness doesn't give one enough reason to neglect some basic concepts concerning with Fourier theory, like the difference between Fourier series and Fourier transform. otherwise, your further deduction probably in the wrong track.

fft: fast Fourier transform for 1-D vector
Note, to really utilize the virtue of fft, the size of vector must a power of 2. otherwise it will fall back to the general discrete Fourier transform
example:
u = [1, 2, 4, 1];
F = fft(u);

the result is as follows:
F = [8.0000 + 0.0000i  -3.0000 - 1.0000i   2.0000 + 0.0000i  -3.0000 + 1.0000i]

fft2: fast Fourier transform for 2-D matrix
identical to fft, if the size of matrix along some dimension is not the power of 2, along that dimension, it automatically falls back into discrete Fourier transform.
example:
S = fspecial('sobel');
F = fft2(S);

the results are as follows:
S = [

     1     2     1
     0     0     0
    -1    -2    -1
];
F = [
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   6.0000 - 3.4641i  -1.5000 - 0.8660i  -0.0000 + 1.7321i
   6.0000 + 3.4641i  -0.0000 - 1.7321i  -1.5000 + 0.8660i
];

fftshift: shift the DC component to the center in order to benefit inspection

example:
FC = fftshift(F);

the result is as follows:
FC = [

  -1.5000 + 0.8660i   6.0000 + 3.4641i  -0.0000 - 1.7321i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
  -0.0000 + 1.7321i   6.0000 - 3.4641i  -1.5000 - 0.8660i
];

ifftshift: inverse of fftshift

following, i will give a comprehensive example to demonstrate the usage of above functions. the purpose of the example is apply low-pass filter in the frequency domain to some image, so as to observe the results:

function freqfilt(img_path)

try
    img = imread(img_path);
catch
    msgbox('cannot open image file, check the file name!');
    return ;
end

h1 = figure(1); imshow(img);
figure(h1); title('Original image');

% Convert into gray-level image if necessary. you've been warned

% FFT transform
img_freq = fft2(img);

% Shift the DC component into center for observing
img_freq = fftshift(img_freq);

% Adjust the frequency value for displaying
img_freq_mag = abs(img_freq);
img_freq_mod = log(img_freq_mag + 1);

% Normalize the matrix since elements are of floating point datatype

f_max = max(max(img_freq_mod));
img_freq_mod = img_freq_mod / f_max;

h2 = figure(2); imshow(img_freq_mod);
figure(h2); title('Spectrum distribution');

% Create a mask for frequency filed filting
[cx, cy] = size(img);

mask = bsxfun(@plus, ((1:cx) - ceil(cx/2)) .^ 2, (transpose(1:cy) - ceil(cy/2)) .^ 2) < 100^2;
h3 = figure(3); imshow(uint8(mask * 255));
figure(h3); title('Mask image');

% Do the filting
img_freq = img_freq .* mask;

% Shift the freqency back to do inverse fft
img_freq = ifftshift(img_freq);

% Note the relation between ifft by utilizing forwarding fft
img_freq = img_freq';
img_ift = fft2(img_freq);
img_ift = img_ift' / numel(img);

% Trucate to real numbers to surpress round-off error
img_ift = real(img_ift);

h4 = figure(4); imshow(uint8(img_ift));
figure(h4); title('Filtered image');

the results are as follows:



imfilter: filtering a given image with specific kernel
example:
img = imfilter(rgb2gray(imread('test.jpg')), fspecial('guassian'));

note, for all kinds of filters, a normalization is recommended before applying to the image. for low pass filter, it's not necessary for adding the original image for displaying. generally for high pass filter, it's tended to request adding the original image before displaying.

freqz2: frequency response for 2-D matrices.
note, this function is mainly for testing kernels, architect your own utility to display the spectrum of some large image.
example:
freqz2(fspecial('laplacian'));

the result is as follows:

No comments:

Post a Comment