function labels = label_subject(P,training_dir)
% FORMAT labels = label_subject(P,training_dir)
%  P            - Name of NIfTI image file
%  training_dir - name of directory containing training data
%  labels       - name of NIfTI label file

%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run Segmentation to obtain tissue probability maps
%
% Method described in these papers: 
% * Ashburner, J., & Friston, K. J. (2005). Unified
%   segmentation. Neuroimage, 26(3), 839-851.
% * Malone, I. B., Leung, K. K., Clegg, S., Barnes, J.,
%   Whitwell, J. L., Ashburner, J., ... & Ridgway, G. R.
%   (2015). Accurate automatic estimation of total
%   intracranial volume: a nuisance variable with less
%   nuisance. Neuroimage, 104, 366-372.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tpm   = fullfile(spm('dir'),'tpm','TPM.nii');
for k=1:size(P,1),
    preproc.channel(k).vols = {deblank(P(k,:))};
    preproc.channel(k).biasreg = 0.001;
    preproc.channel(k).biasfwhm = 60;
    preproc.channel(k).write = [0 0];
end
preproc.tissue(1).tpm = {[tpm ',1']};
preproc.tissue(1).ngaus = 1;
preproc.tissue(1).native = [1 1];
preproc.tissue(1).warped = [0 0];
preproc.tissue(2).tpm = {[tpm ',2']};
preproc.tissue(2).ngaus = 1;
preproc.tissue(2).native = [1 1];
preproc.tissue(2).warped = [0 0];
preproc.tissue(3).tpm = {[tpm ',3']};
preproc.tissue(3).ngaus = 2;
preproc.tissue(3).native = [0 0];
preproc.tissue(3).warped = [0 0];
preproc.tissue(4).tpm = {[tpm ',4']};
preproc.tissue(4).ngaus = 3;
preproc.tissue(4).native = [0 0];
preproc.tissue(4).warped = [0 0];
preproc.tissue(5).tpm = {[tpm ',5']};
preproc.tissue(5).ngaus = 4;
preproc.tissue(5).native = [0 0];
preproc.tissue(5).warped = [0 0];
preproc.tissue(6).tpm = {[tpm ',6']};
preproc.tissue(6).ngaus = 2;
preproc.tissue(6).native = [0 0];
preproc.tissue(6).warped = [0 0];
preproc.warp.mrf = 1;
preproc.warp.cleanup = 1;
preproc.warp.reg = [0 0.001 0.5 0.05 0.2];
preproc.warp.affreg = 'mni';
preproc.warp.fwhm = 0;
preproc.warp.samp = 3;
preproc.warp.write = [0 0];
out1 = spm_preproc_run(preproc);

% Cleanup
delete(out1.param{1});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Refine the registration with the templates in the
% training directory
%
% Method described in these papers:
% * Ashburner, J., & Friston, K. J. (2011). Diffeomorphic
%   registration using geodesic shooting and Gauss–Newton
%   optimisation. NeuroImage, 55(3), 954-967.
% * Ashburner, J., & Friston, K. J. (2009). Computing average
%   shaped tissue probability templates. Neuroimage,
%   45(2), 333-341.
% * Ashburner, J. (2007). A fast diffeomorphic image
%   registration algorithm. Neuroimage, 38(1), 95-113.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
warp1.images{1}(1) = out1.tiss(1).rc;
warp1.images{2}(1) = out1.tiss(2).rc;
warp1.templates = cellstr(spm_select('FPList',training_dir,'^Template_.*\.nii'));
out2 = spm_shoot_warp(warp1);

% Cleanup
delete(out1.tiss(1).rc{1});
delete(out1.tiss(2).rc{1});
delete(out2.vel{1});
delete(out2.jac{1});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform the label fusion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PT     = strvcat(out1.tiss(1).c{1}, out1.tiss(2).c{1});
labels = do_prop(PT,out2.def{1}, training_dir);

% Cleanup
delete(out1.tiss(1).c{1});
delete(out1.tiss(2).c{1});
delete(out2.def{1});

