#!matlab

dataSDE = 

       dwi: [4-D int16]
    tensor: [3x3x131 double]

#!matlab





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% single diffusion encoding (SDE)
% this is a simple clinical scan with 2 shells at b = 1000 and b = 2000
% resolution 1.5mm isotropic

% load the data
x = load('data_examples/single_diffenc');
dataSDE = x.dataSDE;



%% perform estimation without visuliuation
% when baydiff is called for the first time with a specific DWI-scheme 
% it automatically performs the training step (model is saved in /tmp/baydif (linux) ot C:/Temp/baydif (windows) )

[vols model] = baydiff(dataSDE);

% vols is a 4D array and contains the estimated maps, in 
% model.parameter_names the corresponding names are given


%% showing the maps and browsing the signal
% you can directly call browse_micro which also performs estimation 
% and opens immediately a viewer to browse the results.

browse_micro(dataSDE);

% click into the image to show parameters of the voxel.
% use the slider to browse through slices.
% use the popup-menu to the lower left to change contrast
% top/right shows parameters, 
% left to this the spherical average of the log(signal) and the log of the spherical avgerage is shown
% middle/right signal (dots) and prediction (solid lines) are shown (different colors indicate different b-values)
% bottom/right histograms of parameters are shown (diffusivities and volume fractions




%% customize parameters
% via 
% > help baydiff
% you see all the parameters available. The parameters are passed via
% parameter/value pairs. For example, increase the SH order to 4 by 
% > browse_micro(d,d.dwi(:,:,:,1)>150,'noise',20,'lmax',4,'force_retraining',true)
% Note that you have to add the 'force_retraining'-option, because there is
% already a trained model for this scheme





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% multi diffusion encoding
%
% this is data acquired with ordinary SDE and planar MDE encodings (DDE)
% both SDE/DDE are acquired for 6 shells (b=500,1000,1500,2000,2600,3100)
% Data was preprocessed by MPdenoising (Veraat et al) followed by
% deGibbsing (Kellner et al, MRM 2016)


% load the data
x = load('data_examples/multi_diffenc');
dataMDE = x.dataMDE;


%% if your computer is slow/small RAM  do NOT try this
browse_micro(dataMDE,dataMDE.mask, ...
        'force_retraining', true , ...
        'noise',dataMDE.sigma, ...          % noise estimated by MPdenoising
        'traceequality',false, ...          % do not use trace equality constraint for MDE data
        'includeb0',false);                 % do not use b0, if you have enough shells


%% try MDE for a subset of gradients (for small computers)
% only use low high b for estimation

b = squeeze(dataMDE.tensor(1,1,:)+dataMDE.tensor(2,2,:)+dataMDE.tensor(3,3,:));
idx = b<250 | (b>750 & b<1250) | (b>1750 & b<2250);
dataMDE_small = dataMDE;
dataMDE_small.dwi = dataMDE.dwi(:,:,:,idx);
dataMDE_small.tensor = dataMDE.tensor(:,:,idx);
browse_micro(dataMDE_small,dataMDE.mask, ...
        'force_retraining', true , ...
        'noise',dataMDE.sigma, ...          % noise estimated by MPdenoising
        'traceequality',false, ...          % do not use trace equality constraint for MDE data
        'includeb0',true);                 % here we have to use b0, otherwis this is just not enough information





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% fast SDE acquisition
% this is an putative acute stroke protcol (<2min acq time) with an hexagonal 
% q-space b<=2000 at resolution 1x1x5 mm , overall 28 diffusion weightings

% load data from niftis with FSL-style bvecs/bvals
dataStroke = bd_loadData_nii('data_examples/data_hex28.nii.gz',{'data_examples/data_hex28.bval','data_examples/data_hex28.bvec'});

% visualize/browse fitted data
browse_micro(dataStroke,[], ...
                            'qspace','qball', ...      % for non-shell datasets, use qball (radial projections on b^k * e^-b)
                            'force_retraining',true);  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% save estimates to nifti

% estimate
[vols model] = baydiff(dataStroke);

% vols is a 4D array and contains the estimated maps, in 
% model.parameter_names the corresponding names are given

% note that there has to be a nifti present serving as a template
% for saving (in the code below named as dataStroke.nii)

% and now save the stuff

outp = '.';
for k = 1:11,
    tpm = dataStroke.nii;
    tpm.img = single(vols(:,:,:,k));  
    tpm.hdr.dime.datatype = 16;
    tpm.hdr.dime.bitpix = 32;    
    tpm.hdr.dime.dim(5) = 1;
    tpm.hdr.dime.scl_slope = 1;
    tpm.hdr.dime.scl_inter = 0;

    finame = fullfile(outp,[model.parameter_names{k} '.nii'] );
    save_untouch_nii(tpm,finame);
end;

#!matlab

[d.dwi d.noise] = MPdenoising(d.dwi,[],[],'full');

d.dwi = single(ringRm(double(d.dwi),[1 3 10]));


% note that you have to compile the unring code
mex -lfftw3 ringRm.cpp