function Hp = generate_operators(ni,ranks)
% function Hp = generate_operators(ni,ranks)
%
% Generates a family of orthogonal operators
% for patches of size ni x ni.
% The operators are derived by assuming Gaussian
% point spread functions.
%
% 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
%
auto_rank = 0;
if nargin<2
   % if no rank is provided, determine it from PSF
   auto_rank = 1;
end
if nargin<1
   error('Too few arguments!');
end
%%%%%%%%%%%%%%%%%%%%%
% image coordinates %
%%%%%%%%%%%%%%%%%%%%%
[X,Y] = meshgrid([-(ni-1)/2:(ni-1)/2],[-(ni-1)/2:(ni-1)/2]);
x1 = repmat(X(:),1,ni*ni);
x2 = repmat(Y(:),1,ni*ni);
Dx = x1-x1';
Dy = x2-x2';
tolerance = 1e-5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parameters of the optics %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = 35e-3; % focal length
Fnumb = 4; % F-number
D = F/Fnumb; % aperture (actual 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 orthogonal operators %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Computing orthogonal operators (Gaussian kernel)\n');
fprintf('%i Levels\n    ',numlevels)
for k=1:numlevels
   fprintf('\b\b\b\b[%.2i]',k);
   Depth = depthlevel(k);
   for j1=1:length(p)
      % blurring radius
      sigma1 = (gamma*D/2*abs(1-p(j1)*(1/F-1/Depth))).^2;
      % integrate over pixel area
      if (sigma1>tolerance)
          wsi = Dx/sqrt(2)/sigma1;
          wsj = Dy/sqrt(2)/sigma1;
          wdx = .5/sqrt(2)/sigma1;
          derfx = erf(wsi+wdx)-erf(wsi-wdx);
          derfy = erf(wsj+wdx)-erf(wsj-wdx);
          PSF1 = .25*derfx.*derfy;
      else
          PSF1 = (Dx==0).*(Dy==0); % Dirac delta
      end
      % compute a defocused image for each focus setting
      for j2=1:length(p)
         % blurring radius
         sigma2 = (gamma*D/2*abs(1-p(j2)*(1/F-1/Depth))).^2;
         % integrate over pixel area
         if (sigma2>tolerance)
             wsi = Dx/sqrt(2)/sigma2;
             wsj = Dy/sqrt(2)/sigma2;
             wdx = .5/sqrt(2)/sigma2;
             derfx = erf(wsi+wdx)-erf(wsi-wdx);
             derfy = erf(wsj+wdx)-erf(wsj-wdx);
             PSF2 = .25*derfx.*derfy;
         else
             PSF2 = (Dx==0).*(Dy==0); % Dirac delta
         end
         M((j1-1)*ni*ni+[1:ni*ni],(j2-1)*ni*ni+[1:ni*ni]) = ...
             PSF1*PSF2';
      end
   end
   HHT(:,:,k) = M;
end
fprintf('Done.\n');
Hp = zeros(ni*ni*2,ni*ni*2,numlevels,length(ranks));
% extract regularized orthogonal operator
for focus=1:numlevels
   [u,s,v] = svd(HHT(:,:,focus));
   if auto_rank
      ranks = rank(s);
   end
   for rnk=1:length(ranks)
      Hp(:,:,focus,rnk) = ...
          u(:,ranks(rnk):size(u,2))*u(:,ranks(rnk):size(u,2))';
   end
end
% save all data
save Operators Hp
return