function make_template_auto_restart(make_template_auto_restart_params)
load(make_template_auto_restart_params);
if ~exist('workdir','var')
    error('workdir undefined')
elseif ~exist(workdir,'dir')
    error(['INVALID workdir: ',workdir])
else   
    cd(workdir);
end
for i=i:l
    disp(['Working on image ',num2str(i),'/',num2str(l)])
    workdir=pwd;
    save make_template_auto_restart_params;
    if estimate
        if (NRMsourceweight | SEGsourceweight | ApplyMask)
            [pth,nam,ext,patid] = GetPatientID(V(i).fname);
        else
            [pth,nam,ext] = fileparts(V(i).fname);
        end
        % get the patient ID from the filename
        % apply manually edited brain mask to improve the normalisation.
        if NRMsourceweight
            VWF = fullfile(pth,'edits','bin', ['MSK_',patid,'.img']);
            if ApplyMask, msk=spm_vol(VWF); end
        end
        
        % in mls_normalise_DEV_2:
        % VWG       - template weighting image
        % VWF       - source weighting image
        %VWF = dnrm.estimate.weight;
        if  SEGsourceweight
            dseg.estimate.affreg.sourceweight =fullfile(pth,'edits','bin', ['MSK_',patid,'.img']);
            if ApplyMask, msk = spm_vol(dseg.estimate.affreg.sourceweight); end
        end
        matname0 = fullfile(pth,[nam '_sn.mat']);
        if exist(matname0,'file')
            disp(['deleting ',matname0])  
            delete(matname0)
            %            m0=load(matname0);
        end
        
        if optimal
            
            % Estimate spatial normalisation parameters
            % by matching GM with a GM template
            disp(['segmenting ',nam,ext])
            if exist('m0','var') & isfield(m0,'M')
                disp(['using initial M estimate from ',matname0])
                sflags=dseg;
                rerun=1;
                while rerun<4 & rerun
                    %disp(['rerun3 = ',num2str(rerun)])
                    [VT,V(i),rerun,sflags] = spm_segment_mls2(V(i),m0.M,sflags,rerun);
                end
            else
                sflags=dseg;
                rerun=1;
                while rerun<4 & rerun
                    %disp(['rerun4 = ',num2str(rerun)])
                    [VT,V(i),rerun,sflags] = spm_segment_mls2(V(i),VG0,sflags,rerun);
                end
            end
            if optimal==1
                VT(1).fname = fullfile(pth,['tmpGM_'  nam ext]);
                spm_write_vol(VT(1),double(VT(1).dat));
                VTname=VT(1).fname;
            elseif optimal==2
                VT(2).fname = fullfile(pth,['tmpWM_'  nam ext]);
                spm_write_vol(VT(2),double(VT(2).dat));
                VTname=VT(2).fname;
            end    
        else
            VTname=V(i).fname;
            VG1 = VG0;
        end
        clear VT
        VT=spm_vol(VTname);
        % if (ApplyMask)
        %     % apply manually edited brain mask to improve the segmentation
        %     VTo  = ApplyBrainMask(VT,msk);
        %     clear VT
        %     VT = VTo;
        %     clear VTo
        % end
        matname1 = fullfile(pth,[nam '_snc.mat']);
        disp(['normalizing ',VTname])
        if  exist('matname0','var') & exist(matname0,'file')
            prm = mls_normalise_DEV_2(VG1,VT,matname1,VWG,VWF,dnrm.estimate,matname0);
        else
            prm = mls_normalise_DEV_2(VG1,VT,matname0,VWG,VWF,dnrm.estimate);
        end
        
        clear VT
        if exist('VTname') & optimal
            % remove temporary GM image
            mls_cleanup({VTname});
        end
        % Create a spatially normalised version of the original image & brainmask
        disp(['normalizing Anatomical image ...']);
        VN  = spm_write_sn(V(i),prm,dnrm.write);
        VN.fname = fullfile(pth,['w_' nam ext]);
        
        %% Do not write volumes to disk.
        %  spm_write_vol(VN,VN.dat);
        %	if modulation
        %	  disp(['writing modulated normalized Anatomical image ...']);
        %	  spm_write_sn(VN,prm,'modulate');
        %	end
        
        if ApplyMask
            % save the mask image
            disp(['writing normalized Mask image ...']);
            BrainMaskN  = spm_write_sn(msk,prm,dnrm.write);
            BrainMaskN.fname = fullfile(pth,['w_BrainMask' nam ext]);
            foo = BrainMaskN.fname;
            spm_write_vol(BrainMaskN, BrainMaskN.dat);
            clear BrainMaskN;
            BrainMaskN = spm_vol(foo);
        end
        % Segment the spatially normalised image
        disp(['segmenting ',VN.fname])
        sflags=dseg;
        rerun=1;
        while rerun<4 & rerun
            %disp(['rerun2 = ',num2str(rerun)])
            [VT,VN,rerun,sflags]  = spm_segment_mls2(VN,eye(4),sflags,rerun);
        end
        
        
        %	VT(1).fname = fullfile(pth,['tmpGM2_'  nam ext]);
        %	VT(2).fname = fullfile(pth,['tmpWM2_'  nam ext]);
        %	spm_write_vol(VT(1),double(VT(1).dat));
        %	spm_write_vol(VT(2),double(VT(2).dat));
        %	 VT1name=VT(1).fname;
        %	 VT2name=VT(2).fname;
        VT1=VT(1);%spm_vol(VT1name);
        VT2=VT(2);%spm_vol(VT2name);
        VT_3=VT(3);
        if (ApplyMask)
            disp('applying normalized brain mask to improve segmentation')
            VT_1  = ApplyBrainMask(VT1,BrainMaskN);
            VT_2  = ApplyBrainMask(VT2,BrainMaskN);
            mfname{1} = fullfile(pth,['btmpGM2_'  nam ext]);
            mfname{2} = fullfile(pth,['btmpWM2_'  nam ext]);
            VT_1.dat = uint8(get_voxel_values3d(VT_1,VT_1.dim(3),0));
            VT_2.dat = uint8(get_voxel_values3d(VT_2,VT_2.dim(3),0));
            %  VT_3.dat = get_voxel_values3d(VT_3,VT_3.dim(3),0);
        else
            VT_1 = VT(1);
            VT_2 = VT(2);
        end
        VT_1.fname = fullfile(pth,['Go'  nam ext]);
        VT_2.fname = fullfile(pth,['Wo'  nam ext]);
        VT_3.fname = fullfile(pth,['Co'  nam ext]);
        
        % CSF is better estimated using gray and white matter values
        t1 = ((double(VT_1.dat) + double(VT_2.dat) + double(VT_3.dat)) > 0.5*255);
        t2 = (1 - (double(VT_1.dat) + double(VT_2.dat))/255);
        VT_3.dat = uint8(round((t1.*t2)*255));
        
        %  VT_1.dat = uint8(round(double(VT_1.dat)));
        %  VT_2.dat = uint8(round(double(VT_2.dat)));
        
        if size(VT_1.pinfo,1)>2, VT_1.pinfo(3:end)=[]; end
        if size(VT_2.pinfo,1)>2, VT_2.pinfo(3:end)=[]; end
        if size(VT_3.pinfo,1)>2, VT_3.pinfo(3:end)=[]; end
        
        % if modulation,
        %   disp(['writing modulated segmented images ...']);
        %   spm_write_sn(VT_1,prm,'modulate');
        %    spm_write_sn(VT_2,prm,'modulate');
        %    spm_write_sn(VT_3,prm,'modulate');
        % end
        %if write_seg
        %    disp(['writing unmodulated segmented images ...']);
        %    spm_write_vol(VT_1,double(VT_1.dat)/255);
        %    spm_write_vol(VT_2,double(VT_2.dat)/255); 
        %    spm_write_vol(VT_3,double(VT_3.dat)/255); 
        %end        
    else
        VN = V(i);
        VT = [V1(i) V2(i) V3(i)];
    end
    
    % Subsample to a lower resolution, and add the
    % current images to the totals.
    disp('subsampling and updating totals ...')
    for p=1:VG0.dim(3),
        
        M   = VN.mat\VG0.mat*spm_matrix([0 0 p]);
        tn  = spm_slice_vol(VN   ,M,VG0.dim(1:2),1);
        tg  = spm_slice_vol(VT(1),M,VG0.dim(1:2),1);
        tw  = spm_slice_vol(VT(2),M,VG0.dim(1:2),1);
        tc  = spm_slice_vol(VT(3),M,VG0.dim(1:2),1);
        
        msk = find(finite(tn));
        tmp = n(:,p); tmp(msk) = tmp(msk) + tn(msk); n(:,p) = tmp;
        tmp = g(:,p); tmp(msk) = tmp(msk) + tg(msk); g(:,p) = tmp;
        tmp = w(:,p); tmp(msk) = tmp(msk) + tw(msk); w(:,p) = tmp;
        tmp = c(:,p); tmp(msk) = tmp(msk) + tc(msk); c(:,p) = tmp;
        tmp = s(:,p); tmp(msk) = tmp(msk) + 1;       s(:,p) = tmp;
        
    end;
    if exist('VT1name','var')
        mls_cleanup({VT1name,VT2name}); 
    end
    if exist('mfname','var')
        mls_cleanup(mfname); 
    end
    
end;

clear VN
clear VT

% Write out the averages (dividing the totals by the number
% of finite observations).

VN = struct('fname','',...
    'dim',[VG0.dim(1:3) spm_type('uint8')],...
    'mat',VG0.mat,...
    'descrip','');

n=reshape(n./s,VG0.dim(1:3));
g=reshape(g./s,VG0.dim(1:3));
w=reshape(w./s,VG0.dim(1:3));
c=reshape(c./s,VG0.dim(1:3));

VN.fname = 'template.img' ; VN.descrip = 'Customised Template'; spm_write_vol(VN,n);
VN.fname =     'gray.img' ; VN.descrip =   'Grey matter prior'; spm_write_vol(VN,g);
VN.fname =    'white.img' ; VN.descrip =  'White matter prior'; spm_write_vol(VN,w);
VN.fname =      'csf.img' ; VN.descrip =           'CSF prior'; spm_write_vol(VN,c);

clear VN;
VN = struct('fname','',...
    'dim',[VG0.dim(1:3) spm_type('uint8')],...
    'mat',VG0.mat,...
    'descrip','');
%save
VN.fname = 's8template.img'; VN.descrip = 'Customised Template conv 8 mm'; VN = smooth_vol(VN,n,8); spm_write_vol(VN,VN.dat);
stemplate=[pwd,'/',VN.fname];
VN.fname =     's8gray.img'; VN.descrip =   'Grey matter prior conv 8 mm'; VN = smooth_vol(VN,g,8); spm_write_vol(VN,VN.dat);
spriors{1}=[pwd,'/',VN.fname];
VN.fname =    's8white.img'; VN.descrip =  'White matter prior conv 8 mm'; VN = smooth_vol(VN,w,8); spm_write_vol(VN,VN.dat);
spriors{2}=[pwd,'/',VN.fname];
VN.fname =      's8csf.img'; VN.descrip =           'CSF prior conv 8 mm'; VN = smooth_vol(VN,c,8); spm_write_vol(VN,VN.dat);
spriors{3}=[pwd,'/',VN.fname];

if nargout>0
    varargout{1} = stemplate;
end
if nargout>1
    varargout{2} = strvcat(spriors{1},spriors{2},spriors{3});
end
if nargout>2
    varargout{3} = V;
end
try
    owd=pwd;
    cd(workdir)
    delete make_template_auto_restart_params.mat;
    cd(owd)
end
return

% ------------------------------------------------------------------------------------------------------------------------

function V = smooth_vol(V,Y,s)
% smooth volume with gaussian kernel
% based on spm_smooth

% compute parameters for spm_conv_vol
%-----------------------------------------------------------------------
vx = sqrt(sum(V.mat(1:3,1:3).^2));
s  = s./vx;					% voxel anisotropy
s  = max(s,ones(size(s)));			% lower bound on FWHM
s  = s/sqrt(8*log(2));				% FWHM -> Gaussian parameter

x  = round(6*s(1)); x = [-x:x];
y  = round(6*s(2)); y = [-y:y];
z  = round(6*s(3)); z = [-z:z];
x  = exp(-(x).^2/(2*(s(1)).^2)); 
y  = exp(-(y).^2/(2*(s(2)).^2)); 
z  = exp(-(z).^2/(2*(s(3)).^2));
x  = x/sum(x);
y  = y/sum(y);
z  = z/sum(z);

i  = (length(x) - 1)/2;
j  = (length(y) - 1)/2;
k  = (length(z) - 1)/2;

spm_conv_vol(Y,Y,x,y,z,-[i,j,k]);
V.dat = Y;
return


%_______________________________________________________________________
function nits = get_nits(nits)
% Get number of nonlinear iterations
% FORMAT nits = get_nits(nits)

if prod(nits) > 0,
    tmp2 = [1 3 5 8 12 16];
    tmp = find(tmp2 == nits);
    if isempty(tmp) tmp = length(tmp2); end;
    nits = spm_input(['# Nonlinear Iterations?'],'+1','m',...
        ['1  nonlinear iteration |3  nonlinear iterations'...
            '|5  nonlinear iterations|8  nonlinear iterations'...
            '|12 nonlinear iterations|16 nonlinear iterations'],tmp2, tmp);
else, nits = 0; end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function cutoff = get_cutoff(cutoff)
% Get cutoff frequency of DCT bases
% FORMAT cutoff = get_cutoof(cutoff)

tmp2 = [15 20 25 30 35 40 45 50 60 70 80 Inf];
tmp = find(tmp2 == cutoff);
if isempty(tmp) tmp = length(tmp2); end;
cutoff = spm_input('Cutoff spatial normalization','+1','m',...
    ['15mm cutoff|20mm cutoff|25mm cutoff|30mm cutoff|'...
        '35mm cutoff|40mm cutoff|45mm cutoff|50mm cutoff|'...
        '60mm cutoff|70mm cutoff|80mm cutoff|Affine only'], tmp2, tmp);
return;
%_______________________________________________________________________

%_______________________________________________________________________
function reg = get_reg(reg)
% Get amount of regularisation
% FORMAT reg = get_reg(reg)

tmp2 = [100 10 1 0.1 0.01];
tmp = find(tmp2 == reg);
if isempty(tmp) tmp = length(tmp2); end;
reg = spm_input('Nonlinear Regularization','+1','m',...
    ['Extremely Heavy regularization (100)|Heavy regularization (10)|'...
        'Medium regularization (1)|Light regularization (0.1)|'...
        'Very Light regularization (0.01)'], tmp2, tmp);
return;
%_______________________________________________________________________

function Vo = ApplyBrainMask(V,Mask)
% P =spm_get(Inf,'IMAGE',['Select images to Maskify']);
% V   = spm_vol(P);%spm_get(Inf,'*.IMAGE',['Select ',GroupLabel, ' images to normalize & segment']));
f = 'i1.*i2';
for i=1:length(V)
    %V=spm_vol(P(i));
    [pth,nam,ext] = fileparts(V(i).fname);
    Q = ['b' nam];
    
    %-Work out filename for output image
    %------------------------------------------------------------------
    Qdir = spm_str_manip(Q,'Hv');
    Qfil = [spm_str_manip(Q,'stv'),'.img'];
    if ~exist(Qdir,'dir')
        warning('Invalid directory: writing to current directory')
        Qdir = '.';
    end
    Q = spm_get('CPath',Qfil,Qdir);
    Vo = struct(	'fname',	Q,...
        'dim',		[V(i).dim(1:end)],...
        'mat',		V(i).mat,...
        'descrip',	'spm - algebra');
    tmp(1)=V;
    tmp(2)=Mask;
    Vo = spm_imcalc(tmp,Vo,f);
end

%_______________________________________________________________________
function weight = get_weight(weight,type)
% Get cutoff frequency of DCT bases
% FORMAT cutoff = get_cutoof(cutoff)
if (strcmp(type,'NRM template') | strcmp(type,'SEG template'))
    tmp2 = [0 1];
    choices = 'no mask|generic mask';
elseif (strcmp(type,'SEG source') | strcmp(type,'NRM source'))
    tmp2 = [0 1];
    choices = 'no mask|subject mask';
    
else
    error('incorrect type')
end
tmp = find(tmp2 == weight);
if isempty(tmp) tmp = length(tmp2); end;
weight = spm_input([type ' Weighting Mask'],'+1','m',...
    choices, tmp2, tmp);

return;
%_______________________________________________________________________

% 
% %_______________________________________________________________________
% function mask = get_mask(varargin)
% % Get cutoff frequency of DCT bases
% % FORMAT cutoff = get_cutoof(cutoff)
% if (nargin > 0), pth = varargin{1}, end
% if (nargin > 1), patid = varargin{2}, end
% 
%         mask = fullfile(pth,'edits','bin', ['MSK_',patid,'.img']);
% return;
% %_______________________________________________________________________


%_______________________________________________________________________
function ApplyMask = get_ApplyMask(ApplyMask)
% 
% Get cutoff frequency of DCT bases
% FORMAT cutoff = get_cutoof(cutoff)
tmp2 = [0 1];
choices = 'no, don''t apply mask|yes, apply mask';
tmp = find(tmp2 == ApplyMask);
if isempty(tmp) tmp = length(tmp2); end;
ApplyMask = spm_input(['Apply binary mask after segmentation?'],'+1','m',...
    choices, tmp2, tmp);

return;
%___________________________________________________________________
function mls_cleanup(imagename)
for i=1:length(imagename)
    [pth,nam,ext,ver]=fileparts(imagename{i});
    tmp1=fullfile(pth, [nam '.*' ver]);
    disp(['clearing temporary file: ',imagename{i}])
    delete(tmp1);
end
