%-----------------------------------------------------------------------
% function spmTissueSegment(input_folder, input_volume, mask_folder)
% BATCH SPM SEGMENTATION
% input_folder + input_volume <-   path where the volume to segment is located.
%
% function call: matlab -nodisplay -nosplash -nodesktop -nojvm -r "spmNewSegment('./', 'T1', './mask.nii'); quit"
%-----------------------------------------------------------------------

%function spmTissueSegment(input_folder, input_volume)
function [segmTissue]=spmTissueSegment(input_folder, input_volume, output)
    
    %SPM 8/12

    [~, spm_version, ~] = spm_fileparts(spm('dir'));
    if strcmp(spm_version,'spm12')
        TMP_path = [spm('dir'),'/tpm/TPM.nii'];
        p_tools = 'spatial';
        p_preproc = 'preproc';    
    else 
        TMP_path = [spm('dir'),'/toolbox/Seg/TPM.nii'];
        p_tools = 'tools';
        p_preproc = 'preproc8';
    end

	spm_jobman('initcfg');

	matlabbatch{1}.spm.(p_tools).(p_preproc).channel.vols = {[input_folder input_volume '.nii']};
	matlabbatch{1}.spm.(p_tools).(p_preproc).channel.biasreg = 0.0001;
	matlabbatch{1}.spm.(p_tools).(p_preproc).channel.biasfwhm = 60;
	matlabbatch{1}.spm.(p_tools).(p_preproc).channel.write = [0 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(1).tpm = {[TMP_path,',1']};
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(1).ngaus = 2;
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(1).native = [1 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(1).warped = [0 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(2).tpm = {[TMP_path,',2']};
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(2).ngaus = 2;
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(2).native = [1 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(2).warped = [0 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(3).tpm = {[TMP_path,',3']};
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(3).ngaus = 2;
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(3).native = [1 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(3).warped = [0 0];	
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(4).tpm = {[TMP_path,',4']};
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(4).ngaus = 3;
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(4).native = [0 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(4).warped = [0 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(5).tpm = {[TMP_path,',5']};
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(5).ngaus = 4;
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(5).native = [0 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(5).warped = [0 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(6).tpm = {[TMP_path,',6']};
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(6).ngaus = 2;
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(6).native = [0 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).tissue(6).warped = [0 0];
	matlabbatch{1}.spm.(p_tools).(p_preproc).warp.reg = [0 0.001 0.5 0.005 0.2];
	matlabbatch{1}.spm.(p_tools).(p_preproc).warp.affreg = 'mni';
	matlabbatch{1}.spm.(p_tools).(p_preproc).warp.samp = 3;
	matlabbatch{1}.spm.(p_tools).(p_preproc).warp.write = [1 0];

	spm('defaults','pet');
	spm_clf;

	spm_jobman('run',matlabbatch);
    
    
    %Masking the probability maps
    c1 = load_untouch_nii([input_folder, 'c1', input_volume, '.nii']);
    c2 = load_untouch_nii([input_folder, 'c2', input_volume, '.nii']);
    c3 = load_untouch_nii([input_folder, 'c3', input_volume, '.nii']);
    
    c1.hdr.dime.scl_slope=1;
    c2.hdr.dime.scl_slope=1;
    c3.hdr.dime.scl_slope=1;

    mask = ones(size(c1.img));
    mask(c1.img<0.5 & c2.img<0.5 & c3.img<0.5)=0;
    seg = zeros(size(c1.img));
    [~,seg] = max([c1.img(mask==1) c2.img(mask==1) c3.img(mask==1)],[],2);
    
    res=zeros(size(c1.img));
    res(mask==1)=seg;
    segmTissue=zeros(size(c1.img));
    segmTissue(res==3)=1;
    segmTissue(res==1)=2;
    segmTissue(res==2)=3;
    c1.img=segmTissue;
    save_untouch_nii(c1, [input_folder, output]);
    
    disp(['spmNewSegment: segmentation on volume ', input_volume, ' done!']);
end
