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);
    end
end

%% Test Image
% imgname = 'SL'; xtrue = 255 * phantom(256); 
imgname = 'MRI029'; xtrue = double(dlmread('~/mfiles/testImages/TXT/MRI029_256x256.txt'));
% imgname = 'MRIT2NORM057'; xtrue = double( imread('~/mfiles/testImages/BrainWeb/Normal/T2/IMAGES/MRIT2NORM057.tif') );

mx = max(xtrue(:));
mn = min(xtrue(:));
N = size(xtrue);
Nc = N / 2;
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 = N;
params.Nc = Nc;

%% Sampling mask in the Fourier domain and associated fatrices
% % Uniform under-sampling
% dsfact = 2; PSF = zeros( [nx1 nx2] ); PSF(dsfact: dsfact: end, dsfact: dsfact: end) = ones( [nx1 nx2]/dsfact ); sampname = ['PEby' int2str(dsfact)]; 

% % Radial sampling
% nlines = 30;
% PSF = generate_RadialSamplingMask(nx1, nx2, nlines); sampname = ['Rad' int2str(nlines)];

% % Gaussian random sampling along phase-encode in 2D
% % pelocrand = floor( randn(nx1, 1) * Nc(1) ^ 0.85 + Nc(1) ) + 1; pelocrand(pelocrand < 1) = 1; pelocrand(pelocrand > nx1) = nx1;
% perand = zeros(nx1, 1); perand(pelocrand) = 1; 
% PSF = fftshift(repmat(perand, [1 nx2]));
% lSP = length( find( PSF > 0 ) );
% sampname = [ 'RndPE' num2str(lSP / Npix * 100 ) ]; % Add %undersampling to sampname

% % Poisson-Disk Sampling: variable density sampling
% nlayers = 1; rate = [1 1] * sqrt( 100 / 30 ); betapow = 0.8;
% reduction = max( [nx1 nx2] ) ./ round( max( [nx1 nx2] ) * mylinspace( nlayers, 1, nlayers ) / nlayers )
% dens = log10( mylinspace( betapow ^1, betapow ^( max( rate(:) ) / 1.5 ), nlayers ) ) / log10( betapow )
% patt = cell(nlayers, 1); patt{1} = 'rect';
% PSF = generate_PoissonDisk_VDDisjoint_SamplingMask([nx1 nx2], rate, nlayers, reduction, dens, patt);
% PSF = fftshift( PSF );
% sampname = [ 'PDisk' int2str(nx1) ];
% lSP = length( find( PSF > 0 ) );
% sampname = [ sampname num2str( lSP / Npix * 100, '%0.2f' ) ]; % Add %undersampling to sampname
% save( [ sampname, '.mat' ], 'PSF', 'sampname' );
load('/n/calumet/x/sramani/mfiles/WSURE/PDisk256S40.91.mat', 'PSF', 'sampname');

%% Add full-sampling at the center of kspace
PSF = fftshift( PSF );
if( isempty( strfind( sampname, 'Rad' ) ) )
    csfact = 17; csfact2 = (csfact - 1) / 2;
    PSF( Nc(1) - csfact2 + 1 : Nc(1) + csfact2 + 1, Nc(2) - csfact2  + 1 : Nc(2) + csfact2 + 1 ) = 1;
end
PSF = fftshift(PSF);

PSF2 = PSF;

params.PSF = PSF;
params.PSF2 = PSF2;

Ndata = length( find( PSF > 0 ) );
psamp = Ndata / Npix * 100;
sampname = [sampname '_psamp' num2str(psamp, '%0.2f')];
disp(['% sampling = ', num2str(psamp, '%0.2f') ' %']);

%% Fatrices for the data term
Q = 1 / sqrt( Npix ) * Gdft( 'mask', true( [nx1 nx2] ) );
M = Gdiag( PSF );

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

params.Axtrue =  Axtrue;
params.Ndata = Ndata;
params.Q = Q;
params.M = M;
params.A = A;
params.AA = AA;
params.sampname = sampname;
params.AAr = PSF2;

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

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 '_' sampname];
dirsetupmat = '~/mfiles/WSURE/ZMAT/CSMRI/';
dirsetupeps = '~/mfiles/WSURE/ZEPS/CSMRI/';

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

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