% dip06_4.m              gr_11nov21                             
% Filtering in the transform domain
clear all;
return

%% ----------------------- some well-known filters
% check what happens to signals at (pi,pi), (0,pi), (0,pi/2):
%   + - +        + + +        + + +
%   - + -        - - -        0 0 0
%   + - +        + + +        - - -
%                             0 0 0
%

 h1 = [1 1 1; 1 1 1; 1 1 1];
 h2 = [1 2 1; 2 4 2; 1 2 1];
%  h1 = [0 -1 0; -1 4 -1; 0 -1 0];
%  h2 = [-1 -1 -1; -1 8 -1; -1 -1 -1]; 
%  h2 = [1 2 1; 0 0 0; -1 -2 -1];
subplot(121);  mesh(abs(fftshift(fft2(h1,256,256)))); % note: interpolated responses
subplot(122);  mesh(abs(fftshift(fft2(h2,256,256))));

%% ---------- ideal filters
close(2)
f = double(imread('.\img\cameraman.tif'));
F = fftshift(fft2(f));

[u,v] = meshgrid(1:256,1:256);
z = sqrt( (u-129).^2 + (v-129).^2 ); 
for r = 128:-9:1;     % set of 'ideal' disk-shaped filters
    H = double(z<r);         %  lowpass
    % H = 1 - H;  %    or H = double(z>r);  highpass
    F1 = F;  F1(round(z)==r)=1e6;
    figure(1);  imshow(log(abs(F1)),[]);
    G = F .* H;	    % missing zero-padding generates small errors
    g = real(ifft2(ifftshift(G)));
    minmax = [min(g(:)) max(g(:))]
    power = mean(g(:).^2)       % note: image energy is reduced
    figure(2); imshow(g,[]); pause;
    end;
    
%% ----------------------------------- Gaussian filters
for r = 128:-9:1;         
    H = exp(-z.^2/(2*r.^2));      %  lowpass
    % H = 1 - H;                  %  highpass
    figure(1); mesh(H);
    G = F .* H;	    % missing zero-padding generates small errors
    g = real(ifft2(ifftshift(G)));  % (visible periodicity artifact for small r)
    minmax = [min(g(:)) max(g(:))]
    figure(2); imshow(g,[]); pause;
    end;

%% --------------------------------------------  Laplacian filter 
H = -z.^2;         %  
figure(1); mesh(H);   % beware, this Laplacian formulation has huge gain
G = F .* H;
g = real(ifft2(ifftshift(G)));
minmax = [min(g(:)) max(g(:))]
figure(2); imshow(g,[]);

%%
% H = 1+z.^2;         %  subtraction of the Laplacian, 
  H = 1 + .0002 * z.^2;    %  i.e. unsharp masking
figure(1); mesh(H)
G = F .* H;
g = real(ifft2(ifftshift(G)));
minmax = [min(g(:)) max(g(:))]
figure(1); imshow(f,[0 255]);
figure(2); imshow(g,[0 255]);

%% UM using exact and approx. Laplacians: scarcely visible differences
f = double(imread('.\img\cameraman.tif'));
N = 256; % N = 511; % correct padding should be 256+256-1
F = fftshift(fft2(f,N,N));
[u,v] = meshgrid(1:N,1:N);
z2n = ((u-N/2).^2 + (v-N/2).^2)/N.^2; %  normalized and squared distance

H = 1 + 7 * z2n;    %  i.e. unsharp masking
G = F .* H;
g = real(ifft2(ifftshift(G)));
g = g(1:256,1:256);  % cropping if N>256

h1 = [0 -1 0; -1 4 -1; 0 -1 0]/8;
g1 = f + 2 * filter2(h1,f);

h2 = [-1 -1 -1; -1 8 -1; -1 -1 -1]/16;
g2 = f + 1.8 * filter2(h2,f);

minmax = [min(g(:)) max(g(:)); min(g1(:)) max(g1(:)); min(g2(:)) max(g2(:))]

figure(3); imshow([g g1 g2],[0 255]);
figure(4); imshow([g1-g g2-g],[-10 10]);


return;
