% Fig. 4.38
% Obtaining frequency domain filters from spatial domain

clear all;
close all;
f = imread('Fig0438(a)(bld_600by600).tif');
imshow(f,'initialmagnification','fit');

F = fft2(f);
S = fftshift(log(1 + abs(F)));
imshow(S,[],'initialmagnification','fit');

% Generate the spatial filter
h = fspecial('sobel')'
freqz2(h)
axis off

% Obtain the filter in frequency domain
PQ = paddedsize(size(f));
H = freqz2(h, PQ(1), PQ(2));
imshow(abs(H),[], 'initialmagnification','fit');

% Rearrange the data so that the origin is at top left of the frequency
% rectangle
H1 = ifftshift(H);
imshow(abs(H1),[], 'initialmagnification','fit');

% Generate the filtered image in the spatial domain

% Zero padding the image with 0s by default
gs = imfilter(double(f), h);  
imshow(gs,[], 'initialmagnification','fit');

% Generate the filtered image in the frequency domain
gf = dftfilt(f, H1);
figure; imshow(gf,[], 'initialmagnification','fit');

% Show the absolute values for edge detection
figure; imshow(abs(gs),[], 'initialmagnification','fit');
figure; imshow(abs(gf),[], 'initialmagnification','fit');

% Thresholding to binary images
figure; imshow(abs(gs) > 0.2*abs(max(gs(:))),[], 'initialmagnification','fit');
figure; imshow(abs(gf) > 0.2*abs(max(gf(:))),[], 'initialmagnification','fit');

% Check the difference
d = abs(gs - gf);
max(d(:))
min(d(:))
