% Sets up simulated data and other other variables for running WSURE EXPERIMENTS
%
% Sathish Ramani, University of Michigan, 2011

clear;
clc

%% Setup path to include Jeff's toolbox
if(isempty(strfind(path, '/alg/')))
    try
%        addpath('~fessler/l/src/matlab/alg/');
%        setup;
    catch ME
%        currdir = pwd;
%        cd('~/mfiles/');
%        mypath = pathdef;
%        path(mypath);
%        cd(currdir);
        addpath(genpath(~/mfiles/alg'));
    end
end

% imgname = 'Cameraman'; xtrue = double( imread('cameraman.tif') );
%% Test Image
% imgname = 'Lena'; xtrue = double( imread('~/mfiles/testImages/lena256.tiff') );
% imgname = 'Peppers'; xtrue = double( imread('~/mfiles/testImages/Misc_JPortilla/peppers256.png') );
% imgname = 'House'; xtrue = double( imread('~/mfiles/testImages/Misc_JPortilla/house.png') );

mx = max(xtrue(:));
mn = min(xtrue(:));
N = size(xtrue);
nx1 = N(1);
nx2 = N(2);
Npix = nx1 * nx2;

params.imgname = imgname;
params.mn = mn;
params.mx = mx;
params.xtrue = xtrue;
params.Npix = Npix;
params.N = [nx1 nx2];

%% Blur, fatrices, etc

% % Laplacian nonseparable
% nblur = 15;
% nblurby2 = floor(nblur / 2);
% ii1 = -nblurby2 : nblurby2; ii2 = ii1;
% [ii1 ii2] = ndgrid(ii1, ii2);
% h = 1 ./ (1 + ii1 .^ 2 + ii2 .^ 2); blurname = ['Lap' int2str(nblur)]; 

% Gaussian
nblur = 17;
stddev = 2;
dx = 8 * stddev / (nblur - 1);
x = -4 * stddev : dx : 4 * stddev;
x2 = x .^ 2;
h = exp( -0.5 * x2 / stddev ^ 2); blurname = ['Gau' int2str(nblur)]; 

% % Uniform
% nblur = 9;
% h = ones(1, nblur); blurname = ['Uni' int2str(nblur)]; 

if(min(size(h)) == 1) % check if it's a 1D blur
    h = h(:) * h; % 2D blur by outer product
end

psf = h / sum(h(:)); % Normalize blur to sum to 1

tmp = ( fft2( psf, nx1, nx2 ) ); 
PSF2 = abs( tmp ) .^ 2;

A = Gblur_S( true( [nx1 nx2] ), 'psf', psf, 'type', 'fft,same');
AA = Gblur_S( true( [nx1 nx2] ), 'psffft', PSF2, 'type', 'fft,same'); % Fatrix corresponding to (A' * A)

ei = zeros( [nx1 nx2] ); ei(1,1) = 1;

PSF = A * ei;
PSF = real( fft2( PSF ) ); % Assume psf in image domain is odd sized and symmetric 

params.nblur = nblur;
params.psf = psf;
params.PSF = PSF;
params.PSF2 = PSF2;
params.A = A;
params.AA = AA;
params.ei = ei;
params.blurname = blurname;

%% Simulate noisy data
Axtrue = A * xtrue; % imfilter(xtrue, psf, 'replicate', 'same'); % Use replicated boundary rather than periodic for generating data
sigpow = sum( abs( Axtrue(:) - mean( Axtrue(:) ) ) .^ 2 ) / Npix;

params.Axtrue = Axtrue;
params.sigpow = sigpow;

%% Generate random vector for WSURE calculation: this can be done offline and so keep it outside of the algo
r = ones( numel( xtrue ), 1 );
r( floor( numel( xtrue ) / 2 ) + 1 : numel( xtrue ) ) = -1;
rand('state', 0);
for ii = 1:10
    permindex = randperm( numel( xtrue ) );
    r = r( permindex );
end
r = reshape( r, [nx1 nx2] );
params.r = r;

randn('state', 100);
rG = randn([nx1 nx2]); rG = rG - mean(rG(:)); vrG = sum( abs( rG(:) ) .^ 2 ) / Npix; rG = rG / sqrt(vrG);
params.rG = rG;

%% Save setup
fnamesetup = [imgname '_' blurname];
dirsetupmat = '~/mfiles/WSURE/ZMAT/DEBLUR/';
dirsetupeps = '~/mfiles/WSURE/ZEPS/DEBLUR/';

params.fnamesetup = fnamesetup;
params.dirsetupmat = dirsetupmat;
params.dirsetupeps = dirsetupeps;

% save([dirsetupmat 'Z_Deblur_' fnamesetup '_Experiment.mat'], 'params');