function norm_fluid_auto_restart(norm_fluid_auto_restart_params)
load(norm_fluid_auto_restart_params);
cd(workdir);
for i=i:l
    disp(['Working on image ',num2str(i),'/',num2str(l)])
    workdir=pwd;
    save norm_fluid_auto_restart_params;
    %%%
    if (NRMsourceweight | SEGsourceweight | ApplyMask)
        [pth,nam,ext,patid] = GetPatientID(V(i).fname);
    else
        [pth,nam,ext] = fileparts(V(i).fname);
    end
    matname0 = fullfile(pth,[nam '_sn.mat']);
    if exist(matname0,'file')
      m0=load(matname0);
    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:
    % 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
    if exist('matfiles','var') & ~isempty(matfiles)
      matname1=matname0;
      try
          matname0 = matfiles{i,:};%pth,[nam '_sn.mat']);
      catch
          matname0 = matfiles{i};
      end
      
      m0=load(matname0);	
      matname1 = fullfile(pth,[nam '_sn2.mat']);	  
    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
    disp(['normalizing ',VTname])
    if  exist('matname0','var') & exist(matname0,'file')
      prm = mls_normalise_DEV_2(VG1,VT,matname1,VWG,VWF,dnrm.estimate,matname0);
  elseif  exist('matname0','var') & ~exist(matname0,'file')
      prm = mls_normalise_DEV_2(VG1,VT,matname0,VWG,VWF,dnrm.estimate);
    else
      error('matname0 not defined')
    end
    clear VT
    if optimal, 
        % remove temporary GM image
        mls_cleanup({VTname});
    end
    % Create a spatially normalised version of the original image & brainmask
    disp(['writing normalized Anatomical image ...']);
    VN  = spm_write_sn(V(i),prm,dnrm.write);
    VN.fname = fullfile(pth,['w_' nam ext]);
    if writeanat
      spm_write_vol(VN,VN.dat);
    end
    % make average of warped images
    n = n + double(VN.dat);
    
%    if modulation & writeanat
%        disp(['writing modulated normalized Anatomical image ...']);
%        spm_write_sn(VN,prm,'modulate');
%    end
    if writeepi
      disp(['writing normalized EPI / PET / other images ...']);
      if exist('W','var')
	if iscell(W)
	k=size(W,2);
	for j=1:k
	  %	  disp(['Saving Image ', num2str(j), '/',num2str(k)]);
	  % Print progress information in command window
	  str   = sprintf('Saving Image %3d/%-3d ',j,k);
	  fprintf('\r%-40s: %20s',str,' ')       
	  try www=W{i,j};
	  catch www=W{i};
	  end
	  [Wpth,Wnam,Wext] = fileparts(www.fname);
	  WN = spm_write_sn(www,prm,dnrm.write);
	  WN.fname = fullfile(Wpth,['w_' Wnam Wext]);
	  spm_write_vol(WN,WN.dat);
	  delname{j} =  fullfile(Wpth,[Wnam Wext]);
	  clear WN
	end
	elseif isstruct(W)
	  [k1,k2]=size(W);
	  if k1==length(V)
	    k=k2;
	    for j=1:k
	      %    disp(['Saving Image ', num2str(j), '/',num2str(k)]);
	      % Print progress information in command window
	      str   = sprintf('Saving Image %3d/%-3d ',j,k);
	      fprintf('\r%-40s: %20s',str,' ')       

	      [Wpth,Wnam,Wext] = fileparts(W(i,j).fname);
	      WN = spm_write_sn(W(i,j),prm,dnrm.write);
	      WN.fname = fullfile(Wpth,['w_' Wnam Wext]);
	      delname{j} =  fullfile(Wpth,[Wnam Wext]);
	      spm_write_vol(WN,WN.dat);
	      clear WN
	    end
	  elseif k2==length(V)
	    k=k1;
	    for j=1:k
	      %     disp(['Saving Image ', num2str(j), '/',num2str(k)]);
	      % Print progress information in command window
	      str   = sprintf('Saving Image %3d/%-3d ',j,k);
	      fprintf('\r%-40s: %20s',str,' ')       
	      
	      [Wpth,Wnam,Wext] = fileparts(W(j,i).fname);
	      WN = spm_write_sn(W(j,i),prm,dnrm.write);
	      WN.fname = fullfile(Wpth,['w_' Wnam Wext]);
	      spm_write_vol(WN,WN.dat);
	      delname{j} =  fullfile(Wpth,[Wnam Wext]);
	      clear WN
	    end
	  end
	  mls_cleanup(delname);
	  clear delname
	end
      end
      fprintf('%s%20s',repmat(sprintf('\b'),1,20),'... done',sprintf('\n'))%-#
    end
    if ApplyMask
      % save the mask image
      disp(['writing normalized Mask image ...']);
      BrainMaskN  = spm_write_sn(msk,prm,dnrm.write);
      BrainMaskN.fname = fullfile(pth,['wBrainMask' 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=spm_vol(VT1name);
    VT2=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,['Gw_'  nam ext]);
    VT_2.fname = fullfile(pth,['Ww_'  nam ext]);
    VT_3.fname = fullfile(pth,['Cw_'  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));
    
    % calculate volumes of GM, WM and CSF
    factor = abs(det(VT_1.mat(1:3,1:3)))/1000;
    gray  = factor*sum(VT_1.dat(:))/255;
    white = factor*sum(VT_2.dat(:))/255;
    csf   = factor*sum(VT_3.dat(:))/255;
    
    % save volumes in log-file
    fprintf(fid,'\n%s\t%5.3f\t%5.3f\t%5.3f\t%5.3f',nam,gray,white,csf);
    
    %  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 writeseg
      if modulation,
        disp(['writing modulated segmented images ...']);
        if max_seg>0, spm_write_sn(VT_1,prm,'modulate');end
        if max_seg>1, spm_write_sn(VT_2,prm,'modulate');end
        if max_seg>2, spm_write_sn(VT_3,prm,'modulate');end
      end
      disp(['writing unmodulated segmented images ...']);
      if max_seg>0,spm_write_vol(VT_1,double(VT_1.dat)/255);end
      if max_seg>1,spm_write_vol(VT_2,double(VT_2.dat)/255);end
      if max_seg>2,spm_write_vol(VT_3,double(VT_3.dat)/255);end
    end
    mls_cleanup({VT1name,VT2name}); 
    if exist('mfname','var')
        mls_cleanup(mfname); 
    end
    clear VT
end    
%%%%

fclose(fid);

% save average image
n = n/length(V);
VN.fname   = 'avgT1.img';
VN.descrip = 'Average of warped T1 images';
spm_write_vol(VN,n);
try
    delete norm_fluid_auto_restart_params.mat;
end
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
return
%_______________________________________________________________________
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;
