% This program helps in the empirical estimation of the optimal
% regularization parameter. Typically, a regularization parameter of 0.1
% works fine. However, if that does not seem satisfactory, then this program
% will help in selecting a suitable regularization parameter.


% Which Input signal would you like?
% (1) Cameraman (default)
% (2) Lenna 
% (3) Theater
% (4) Boats
% (5) Birthday
% The images are grayscale images
Input_signal = 1;

% Which blurring function would you like?
% (1)  Box car blur (9x9 blur) (default)
% (2)  Adjustable Box car blur (you specify the amount of blur size)
% (3)  Circular blur with radius 7
% (4)  Low pass filter used by Kalifa et al. in 
%      "Image Deconvolution in Mirror Wavelet bases", ICIP 98, pg. 565-569,
%      Chicago, 1998.

Blur = 1;

% What blurred-signal-to-noise ratio (in dB) would you like?
% If y = x*h + n is the deconvolution model, the BSNR = |norm(x*h)/norm(n)|^2
% default is 40 dB (* denotes convolution)

BSNR = 40;

% Two wavelet basis are required for implementing the wavelet-domain wiener
% filter (WDWF) used to perform the estimation. 

Wavelet_Filt_type1 = daubcqf(6,'min'); 
Wavelet_Filt_type2 = daubcqf(2,'min'); 

% What threshold should be used in for hard thresholding (1st stage in WDWF)?
% default threshold is 3*sigma, sigma is standard deviation of noise.

Thresh_Factor = 3.0;

% How many levels of decomposition do you want?
% The maximum possible is log(N), however 4 levels are often sufficient.

Decomp_levels = 3;

% Do you want to estimate the power spectral density (PSD) of the input
% signal from the observation or do you want to use the original signal to
% get the PSD? The PSD is required for the Wiener filter and the WaRD
% algorithm.
% Choices are PSD = 'Estimate'
%             PSD = 'Original'
% NOTE: For arbitrary images, the estimation 

PSD = 'Estimate';

% Do you want to estimate the noise variance of the input
% signal from the observation or do you want to use the original signal to
% get the PSD? The noise variance is required for the Wiener filter and the 
% WaRD algorithm.
% Choices are Variance = 'Estimate'
%             Variance = 'Original'
%
% NOTE: PSD estimation is not very robust and might require some tweaking
% for arbitrary images.

Variance = 'Estimate';
Reg_Param = [0.00001 0.0001 0.001 0.01 0.025 0.05 0.07 0.1 0.125 ...
      0.15 0.2 0.25 0.35 0.45 0.6 0.9];

WaRD_Setup
for reg_counter = 1:length(Reg_Param)
  reg_param = Reg_Param(reg_counter);  
  Run_WaRD
  Result_norm(reg_counter) = norm(x_ward(:));
  disp(strcat('Counter = ', num2str(reg_counter),'of',  num2str(length(Reg_Param))))
  Actual_mse(reg_counter) = norm(x_ward(:) - X(:));
end  
figure(2) 
%subplot(2,1,1)
plot(Result_norm,'+')
xlabel('Regularization parameter index')
ylabel('Norm of the WaRD estimate')
title('Choose reg parameter corresponding to the start of the flat part of curve')
%subplot(2,1,2)
%plot(Actual_mse,'+')
%xlabel('Regularization parameter')
%ylabel('Actual MSE of WaRD estimate')

disp('From figure 2, choose the regulazation parameter index from the x-axis where the flat part of the curve begins')
opt_num = input('Enter the reg parameter index. This should be a +ve integer    between 1 and 14\n');
disp(strcat('Near-optimal regularization parameter is =',num2str(Reg_Param(opt_num)),'. Plug this regularization parameter in the Expt_Setup.m program and run the WaRD algorithm.'))

