function Hp = learn_operators(ni,ranks,training_seq)
% function Hp = learn_operators(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 %
%%%%%%%%%%%%%%%%%%%%%
X = ones(ni,1)*[-(ni-1)/2:(ni-1)/2];
Y = [-(ni-1)/2:(ni-1)/2]'*ones(1,ni);
%%%%%%%%%%%%%%%%%%%%%
% 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;
pad = 3;
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)
      TS = [];
      DepthMap = depthlevel(j);
      for i=1:size(training_seq,2)
         Rad = abs(randn(ni+2*pad))*200+100;
         J = zeros(ni+2*pad,ni+2*pad,length(p));
         % use heat equation to compute convolution (Gaussian PSF)
         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
         J = J(pad+[1:ni],pad+[1:ni],:);
         TS = [TS J(:)];
      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