% dip06_1.m                gr16mar15   22oct18
% 2-D DFT calculation and examples
clear all;

%%%%%%% note
iptsetpref('ImshowBorder','tight');
iptsetpref('ImshowInitialMagnification','fit');
%%%%%%%%

N = 16;          %%% reminder: matrix calculation of 1-D DFT
f = rand(N,1);

A = dftmtx(N);   % build the DFT matrix
Ainv = conj(A) ./ N;    % and its inverse

clf
for k=1:N/2+1
    subplot(N/2+1,2,1+2*(k-1)); stem(0:N-1, real(A(k,:)));
    subplot(N/2+1,2,2+2*(k-1)); stem(0:N-1, imag(A(k,:)));
end

F = A * f;      % direct transform
f_hat = Ainv * F;   % inverse transform
err = mean((f-f_hat).^2)

return

%% "separable" and matrix calculations of 2-D DFT
clear all;

M = 30; N = 40;
f = rand(M,N);      % try several times, errors change slightly

Ar = dftmtx(M);  Ac = dftmtx(N);   % build the DFT matrices
Ar_inv = conj(Ar)./M;  Ac_inv = conj(Ac)./N;   % and their inverses

Fc = [];    %%%% "separable" transform
for i = 1:N,  % column-wise transf.
    Fc = [Fc  Ar*f(:,i)];
end;
F = [];
for i = 1:M,  % row-wise transf.
    F = [F;  Fc(i,:)*Ac];
end;

f_hat = ifft2(F);  % standard Matlab inverse transform
err_sep = mean(mean((f-f_hat).^2))

F_matr = Ar * f * Ac;   % matrix formulation of direct...
f_matr = Ar_inv * F_matr * Ac_inv;  % ...and inverse transform
err_matr = mean(mean((f-f_matr).^2))


%% some examples
clear all; clf;

L = 100;  %%%  reminder: DFT of a 1-D rectangular signal
for k=1:round(L/11):L
    x = [ones(1,k)*.9  zeros(1,L-k)];
    X = abs(fft(x));
    subplot(211); plot(x);
    subplot(212); plot(X); axis 'tight'; pause();
    end
return;

%%
clf;
f = zeros(256,256);   f(118:138,108:148)=255;  %  rectangle
imshow(f); pause;
Fm = abs(fft2(f));      % 2-D fft
imshow(log(Fm+10),[]); pause;  % log transform for better visibility
Fm = fftshift(Fm);      % swap opposite quadrants
imshow(log(Fm+10),[]); pause;    
plot(Fm(:,129));  %   1-D plots
    % peak ampl. is M * N * white * relative size of white area
max(Fm(:))
256 * 256 * 255 * (21*41)/(256*256)

%%
f = zeros(256,256);   f(128,128)=255;  % impulse
imshow(f,[]);  pause;
Fm = fftshift(abs(fft2(f)));
imshow(Fm);
%%
f = zeros(256,256);   f(128,:)=255;  % vert. impulse
imshow(f,[]);  pause;
Fm = fftshift(abs(fft2(f)));
imshow(Fm,[]);
%%
f = f-mean(f(:));               % zero-mean impulse
Fm = fftshift(abs(fft2(f)));
imshow(Fm,[]);
%%
f = diag(zeros(1,256)+255);     % diagonal impulse
imshow(f,[]); pause;
Fm = fftshift(abs(fft2(f)));
imshow(Fm,[]);
%%
f = [zeros(256,56)  zeros(256,200)+255];   % horiz. step
imshow(f,[]);  pause;
Fm = log(fftshift(abs(fft2(f)))+.01);
imshow(Fm,[]); pause;
plot(Fm(129,:));   %   1-D plot
%%
f = zeros(256,256)+255;  f = triu(f);       % diag. step
imshow(f,[]);  pause;
Fm = fftshift(abs(fft2(f)));
imshow(Fm.^0.3,[]);              % power map for better visibility
%%
u = .4*pi;  v = .1*pi;
[x,y] = meshgrid([0:1:100],[0:1:100]);
f = real(exp(j*(u*x+v*y)));      % cosine grating
imshow(f,[]);   pause;
Fm = fftshift(abs(fft2(f)));
imshow(Fm,[]);
%%
u = .07*pi;  v = .03*pi;
[x,y] = meshgrid([0:1:100],[0:1:100]);
f = round(real(exp(j*(u*x+v*y)))/2+0.5)-0.5;      % square wave grating
imshow(f,[]);   pause;
Fm = fftshift(abs(fft2(f)));
imshow(Fm,[]); pause();
imshow(Fm.^0.3,[]); % small terms due to periodic border effect are visible
%%
 f = double(imread('.\img\lenar.tif'));     % real image
% f = double(imread('.\img\cameraman.tif'));
figure(1); imshow(f,[]);
Fm = log(fftshift(abs(fft2(f)))+.01);
figure(2); imshow(Fm,[]); 
%%
mask = ones(256,1) * [linspace(0,1,30) ones(1,196) linspace(1,0,30)];
mask = (mask + mask')/2;
 f = f .* mask;                      %  reduce image border artifacts
figure(1); imshow(f,[]);
Fm = log(fftshift(abs(fft2(f)))+.01);
figure(2); imshow(Fm,[]); 

%%
 f = double(imread('.\img\cameraman.tif'));
F = fftshift(fft2(f));                      %  direct and inverse transforms
g = real(ifft2(ifftshift(F)));          % better use 'real' against roundoff errors
err = f-g;
imshow(err,[]);
max(abs(err(:)))   

%%
f = double(imread('.\img\cameraman.tif'));
Fp = fftshift(angle(fft2(f)));         % phase of the DFT
max(Fp(:)), min(Fp(:))
imshow(Fp,[]);  

%%
Fm = abs(fft2(f));  Fp = angle(fft2(f));    %  DFT and ...
 fr = real(ifft2(Fm.*exp(j*Fp)));           %    ... IDFT
% fr = real(ifft2(exp(j*Fp)));              % IDFT, phase only (unitary magnitude)

%  Fp1 = ones(256,1) * linspace(-pi,pi,256);
%  Fp1 = (Fp1 + Fp1')/2;               % IDFT, magnitude only (linear phase)
%  Fp1 = randn(256,256)*4; % +Fp;      % IDFT, magnitude only (random phase)
%  mesh(Fp1); pause;
%  fr = real(ifft2(Fm.*exp(j*Fp1)));
 imshow(fr,[]);

return;
