function Hp = learn_operators_overlap(ni,ranks,training_seq)
% function Hp = learn_operators_overlap(ni,ranks,training_seq)
%
% Learns a family of orthogonal operators
% for patches of size ni x ni from defocused images.
%
% Copyright 2006 Paolo Favaro (p.favaro@hw.ac.uk)
% 
% School of Engineering and Physical Sciences
% Heriot-Watt University, Edinburgh, UK
% 
% Last revision: May 2006
%
% This program can be used only for research purposes.
% This program is distributed WITHOUT ANY WARRANTY; 
% without even the implied warranty of MERCHANTABILITY 
% or FITNESS FOR A PARTICULAR PURPOSE.

%
% ni = square patch size
%
if nargin<3
   need_training = 1;
   training_seq = zeros(ni*ni*2,ni*ni*4);
else
   need_training = 0;
   error('Not yet implemented');
end
auto_rank = 0;
if nargin<2
   % if no rank is defined, then compute rank of training set
   auto_rank = 1;
end
if nargin<1   
   error('Too few arguments!');
end
%%%%%%%%%%%%%%%%%%%%%
% image coordinates %
%%%%%%%%%%%%%%%%%%%%%
m = ni*3;
n = ni*ni*5;
X = ones(m,1)*[-(n-1)/2:(n-1)/2];
Y = [-(m-1)/2:(m-1)/2]'*ones(1,n);
%%%%%%%%%%%%%%%%%%%%%
% optics parameters %
%%%%%%%%%%%%%%%%%%%%%
F = 35e-3; % focal length
Fnumb = 4; % F-number
D = F/Fnumb; % aperture (effective lens diameter)
gamma = .8e4; %calibration parameter (CCD pixel size)
z0 = .53; % distance of near focal plane from camera
z1 = .85; % distance of far focal plane from camera
p = 1./(1/F-1./[z0 z1]);% distance of lens from CCD (two images)
%%%%%%%%%%%%%%%%%%%%%%%
% synthetic depth map %
%%%%%%%%%%%%%%%%%%%%%%%
numlevels = 50; % number of operators (i.e., number of depth
                 % levels we distinguish
% positions of equifocal planes
depthlevel = [z0:(z1-z0)/(numlevels+1):z1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate defocused images %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
randn('state',0);
Nit = 5;
if need_training
   fprintf('Computing orthogonal operators (learning)\n');
   fprintf('%i Levels\n    ',numlevels)
   for j=1:numlevels
      fprintf('\b\b\b\b[%.2i]',j);
      % generate training set (Gaussian PSF)
      DepthMap = depthlevel(j);
      Rad = abs(randn(m,n))*200+100;
      % use heat equation to compute convolution (Gaussian PSF)
      J = zeros(m,n,length(p));
      for k=1:length(p)
          Mask = (gamma*D/2*abs(1-p(k).*(1/F-1/DepthMap))).^2;
          J(:,:,k) = del2(defocusUnif(Rad,Mask));
      end
      TS = zeros(ni*ni*2,size(J,2)-2*ni);
      for i=1:size(J,2)-2*ni
         TS(:,i) = reshape(J(ni+[1:ni],ni+[1:ni]+(i-1),:),...
             ni*ni*length(p),1);
      end
      [u,s,v] = svd(TS);
      if auto_rank
         ranks = rank(s);
      end
      Hp(:,:,j) = u(:,2*ni*ni-ranks:ni*ni*2)*...
          u(:,2*ni*ni-ranks:ni*ni*2)';
   end
end
fprintf('Done.\n');
% save all data
save Operators Hp
return