% This function implements MacCallum's Simulated annealing 
% blind deconvolution algorithm
% 
% g   : the blurred image = convolution (f,h) + noise
% f,h : initial guess for f and h
% Nc  : maximum number of cycle iterations
% Ns  : maximum number of scan iterations per cycle
% 
% reference: B. C. McCallum,"Blind Deconvolution by simulated annealing," Opt. Commun.,
%            vol. 75, no 2, pp. 101-105, Feb. 1990
%
% Alper Corlu 12/11/2002
% corlua@dept.physics.upenn.edu

function [f,h, Qrecord, Trecord] = simanneal_main (g, f, h, Nc, Ns)

% start timing
tic

% reset random number generator state
  rand('state',sum(100*clock))

% STEP 1
% Scale g so that Energy(g) = Energy (conv2(fo,go))
  
  scale_param = sqrt ( energy(g) / energy (conv2(f,h)));
  g = g ./scale_param;

% initialize the counters
  nc = 1; % initialize cycle counter 
  ns = 1; % initialize scan counter

% initialize temperature  
  T = cost_simannealing (f,h,g) / 10;
  T = T / 0.8; % to get above T in step 2 
  fprintf (1,'\nProgress [');

while nc <= Nc
   ns = 1;
   
   % record the values of Q and T at each cycle
     Qrecord(nc) = cost_simannealing (f,h,g);
     Trecord(nc) = T;
     
   % progress bar
     fprintf (1,'.');
   
   % STEP 2
   % Calculate values for T (temperature) and alpha (pertubation scales)
     T = T * 0.8;
     alpha = 100 * sqrt(T);

   % STEP 3
   % Scale the estimates of f and h so that scaled images have equal rms
     [f,h] = scale_rms (f,h);
   
   % STEP 4-5
   % Start the Scans
   
   while ns <= Ns
      
      % first try to perturb f
        f = perturb (f, h, g, alpha, T);
      
      % then for h
        h = perturb (h, f, g, alpha, T);
      
    
      % STEP 6
      % increment scan counter and retrun to STEP 4
        ns = ns +1;
   end
   	
   % STEP 7
   % increment cycle counter and return to STEP 2
     nc = nc +1;
end

fprintf (1,']\n');

% stop timing 
toc