% dip06_5.m           gr18mar15     24oct17     28oct19
% Separable filters and filter design
clear all;
return
%% ------------------------ building separable filters

 b1 = [1 0 -1];
 b2 = [1 2 1];  % vert. Sobel mask

%  b1 = [1 2 3 5 8 5 3 2 1]; 
%   b2 = [1 2 3 5 8 5 3 2 1];
%   b2 = [1 1 1 1 1 1 1 1 1]; 
%   b2 = [-1 -1 -1 -1 8 -1 -1 -1 -1]; % note ripple in freq response
%  b2 = [1 -2 3 -6 8 -6 3 -2 1];  %  ripple in imp. resp --> less ripple in freq
%  b2 = [0 0 0 -1 2 -1 0 0 0];
%  b2 = [0 0 0 -1 3 -1 0 0 0];

h = b1' * b2
DCgain = sum(h(:))
% H = abs(fftshift(fft2(h,256,256)));  mesh(H);
 figure(1); freqz2(h);

f = double(imread('.\img\cameraman.tif'));
g = filter2(h,f);
 figure(2); imshow(g,[]); 


%% check whether a given filter h is separable
h = [1 2 1; 2 4 2; 1 2 1];    % try changing e.g. the center coefficient...

%  h = [1 2 1; 0 0 0; -1 -2 -1];  % warning: ill-conditioned

% "random" 5x5 filter having half-quadrant symmetry
%  z = randn(10,1); z = filter([1 1 1]/3,1,z);
%  a=z(2); b=z(3); c=z(4); d=z(5); e=z(6); f=z(7);
%  h = [a b c b a; ...
%       b d e d b; ...
%       c e f e c; ...
%       b d e d b; ...
%       a b c b a]; 

disp('==== Separable filter that approximates a given freq. response')
disp('rank(h)=1 if h is separable'), disp(rank(h));
[V,D] = eig(h);  % V = eigenvectors matrix;  D = eigenvalues (diag) matrix
eigenvalues = diag(D)'
Vm1 = inv(V);
disp('this must be h again:'),  disp(V*D*Vm1);    
   % extract the separable terms, taking the eigenvector which 
   % corresponds to the single nonzero eigenvalue.
   % Or, build an approximated separable version of h
[dum,ind] = max(abs(diag(D))); 
b1 = V(:,ind)' * sqrt(D(ind,ind));
b2 = Vm1(ind,:) * sqrt(D(ind,ind));
disp('1-D components:'), b1, b2
hs = b1' * b2;
disp('separable filter similar/equal to h:'),  disp(hs);
disp('check that rank(hs)=1:'), disp(rank(hs));
 figure(1); freqz2(h); % desired freq. response
 figure(2); freqz2(hs);% obtained freq. response
 

%% ---------- 2-D filter design using windowing (fwind1 and fwind2 commands)
Hd = zeros(17,17);
Hd(5:13,5:13) = 1;
Hd(8:10,8:10) = 0;    % bandpass, DC is in (9,9)
figure(1); subplot(131); mesh(Hd);

hr = fwind1(Hd,ones(17,1));  % design using rectangular window 
hb = fwind1(Hd,blackman(17));  % (use fwind2 to use 2-D windows) 
subplot(132); freqz2(hr,[34 34]); axis([-1 1 -1 1  0 1]);
subplot(133); freqz2(hb,[34 34]); axis([-1 1 -1 1  0 1]);

figure(2); subplot(121); mesh(hr); subplot(122); mesh(hb);
DCgain_rectangular = sum(hr(:))
DCgain_blackman = sum(hb(:))    % check if gain 0 at DC is in fact achieved...


%%  ---------- similar, with more isotropic Hd
[f1,f2] = freqspace(19,'meshgrid');  % -1 ... 1 (*pi) in 21 steps
Hd = ones(19);
r = sqrt(f1.^2 + f2.^2);
Hd((r<0.2) | (r>0.5)) = 0; 
figure(1); subplot(131); mesh(Hd);
hr = fwind1(Hd,ones(19,1));
hh = fwind1(Hd,hamming(19));
subplot(132); freqz2(hr); axis([-1 1 -1 1  0 1]);
subplot(133); freqz2(hh); axis([-1 1 -1 1  0 1]);
figure(2); subplot(121); mesh(hr); subplot(122); mesh(hb);
DCgain_rectangular = sum(hr(:))
DCgain_hamming = sum(hh(:))    % check if gain 0 at DC is in fact achieved...


%% -----------------  2-D filter design with the McClellan transform
 b = [1 1 1];  b = b / sum(b);  % note differences on diag. edges
% b = firpm(30,[0 0.05 0.15 0.55 0.65 1],[0 0 1 1 0 0]);  % Parks-McClellan 1D design
h1 = ftrans2(b);
h2 = b' * b;   % compare with the separable filter
figure(1); subplot(121); freqz2(h1);
           subplot(122); freqz2(h2);
f = double(imread('.\img\cameraman.tif'));
g1 = filter2(h1,f);
g2 = filter2(h2,f);
figure(2); imshow([g1 g2],[]);
diff = g1 - g2;
minmaxdiff = [min(diff(:)),  max(diff(:))]    % small differences in many cases...
figure(3); imshow(diff,[]);



