function result = blur_mse (x, y_blurred, y, s_0, noise_variance, blur_type, theta, sigma_1, sigma_2)
% USAGE    : result = blur_mse (x, y_blurred, y, s_0, noise_variance, blur_type, theta, sigma_1, sigma_2)
% FUNCTION : Computes the blur-MSE (oracle of the blur-SURE) of an image y, for a series of values 
%            of the PSF scaling factor. 
% INPUT    : * 'x'. The original (ground-truth) image.
%            * 'y_blurred'. The noiseless, but blurred image.
%            * 'y'. The noisy, blurred image.
%            * 'blur_type'. The PSF 'probed' by the blur-MSE may assume (here) only three 
%              different types:
%               1. isotropic Gaussian
%               2. jinc
%               3. Asymmetric, oriented Gaussian, in which case the orientation angle, theta,
%                  and the length of the two axes, sigma_1, sigma_2, must be provided.
%            * 'noise_variance'.
%            * 's_0'. The center of an interval [s_0-0.5,s_0+0.5] where the unknown scaling factor 
%              is supposed to lie.
% OUTPUT   : 'result'. A 3-column array containing
%               1st column. Values of the PSF scaling factor in the interval [s_0-0.5,s_0+0.5].
%               2nd column. Values of the blur-MSE optimal Tikhonov regularization for that scaling factor.
%               3rd column. Values of the optimal blur-MSE for that scaling factor.
% The best scaling factor estimate corresponds to the index for which result(:,3) is minimum.
%
% DATE     : June 2013
% AUTHOR   : Feng Xue, Thierry Blu, the Chinese University of Hong kong, Shatin, Hong Kong
%            mailto:thierry.blu@m4x.org

wait_bar = waitbar(0, 'Blur-MSE computation. Please wait...');
waitbar(1/4, wait_bar);

x = double(x); y = double(y);  
[M,N] = size(x); X = fft2(x); Y = fft2(y); Y_blurred = fft2(y_blurred);

%% %%%%% key parameter --- regularization %%%%%%%%%%%%%

% first stage
S = design_regularizer (M,N, 2); 
% if (s_0>3.0)
%     S = design_regularizer(M,N, 2.5);
% end

lambda_n = noise_variance*logspace(-5, -1, 50);

s_n = linspace( s_0-0.5, s_0+0.5, 5); % s_n = [1 1.5 2 2.5 3];
s_n = s_n'; sure = s_n;  mse = sure;

waitbar(2/4, wait_bar);
fid = fopen('lambda_s.txt', 'wt');

for i = 1:length(s_n)
    s = s_n(i);
    
    switch blur_type
        case 1
            h_n = gaussian_kernel (255, 255, s) ;
        case 2
            h_n = jinc_function (s) ;
        case 3
            h_n = non_symmetric_gaussian_kernel (255, 255, theta,  sigma_1, sigma_2, s ) ;
    end
    
    h_n = centering_kernel (h_n, M, N);      H_n = real(fft2(h_n));      HH_n = H_n.*H_n;

for j = 1:length(lambda_n)
    lambda = lambda_n(j);      U= HH_n./(HH_n+lambda.*S);   
    U1 = U-1; 

% evaluate blur-MSE
dif1 = U.*Y-Y_blurred; mse(j) = real(dif1(:)'*dif1(:))/(M*N*M*N);  

end

lambda_opt_oracle = lambda_n(find(mse==min(mse)));
yy = [s, lambda_opt_oracle(1)];
    
fprintf(fid, '%2.4f %1.4e \n', yy);

end

fclose(fid);

waitbar(3/4, wait_bar);

%% second stage

clear mse s_n lambda_n;

load lambda_s.txt;    aa = lambda_s;   
aa_1 = aa(:,1);   aa_2 = -log10(aa(:,2));   
s_n = linspace(s_0-0.5, s_0+0.5, 50);   s_n = s_n';  mse = s_n;
lambda_n = interp1 (aa_1, aa_2, s_n);  lambda_n = 10.^(-lambda_n);  

for i = 1:length(s_n)
    lambda = lambda_n(i);    s = s_n(i);
    
    switch blur_type
        case 1
            h_n = gaussian_kernel (255, 255, s) ;
        case 2
            h_n = jinc_function (s) ;
        case 3
            h_n = non_symmetric_gaussian_kernel (255, 255, theta,  sigma_1, sigma_2, s ) ;
    end
        
h_n = centering_kernel (h_n, M,N);      H_n = real(fft2(h_n));      HH_n = H_n.*H_n;
U = HH_n./(HH_n+lambda.*S);   U1 = U-1;

% evaluate blur MSE
dif1 = U.*Y-Y_blurred; mse(i) = real(dif1(:)'*dif1(:))/(M*N*M*N);

end

result = [s_n lambda_n mse];

waitbar(4/4, wait_bar);
close (wait_bar);

return;


