function [template,priors,imgVols]=make_template(parr,modality,ImgVols,optimal)
% make_template.m: template creation for T1 or T2 images with SPM2.
% 
% Please note, with some setup work,
% this can be configured to run over Sun Grid Engine.  (SGE)
% By default SGE mode is disabled, but you can
% un-comment line 12 in vbm.m and line 19 in make_template.m 
% to allow SGE mode to be used. you will need to make some other changes in
% the code also.  Email me with any questions about SGE mode.
%
% SEE ALSO: vbm_ui.m
% Jan. 10, 2006
% Matthew L. Senjem
% Mayo Clinic 
% Senjem.Matthew1@mayo.edu
%
% When using this code, please reference the following paper:
%
% Senjem ML, Gunter JL, Shiung MM, Petersen RC, Jack CR Jr., "Comparison of different methodological implementations of voxel-based
% morphometry in neurodegenerative disease." Neuroimage. 2005 Jun;26(2):600-8. 
sge=0;
spm_defaults;
global defaults;

if ~exist('parr','var')
    sge=0;
    % Uncomment the next line to allow SGE mode to be used.
    % sge = spm_input('Run over SGE?','+1','m','YES|NO',[1 0], 1);
elseif exist('parr','var') && strcmpi(parr,'sge')
    sge=1;
end
if ~exist('modality','var')
    moda = spm_input('Modality?','+1','m','T1|T2',[1 2], 1);
    switch moda
        case 1
            modality='t1';
        case 2
            modality='t2';
    end
end
if ~exist('ImgVols','var')
    ImgVols = spm_vol(spm_get(Inf,'*.IMAGE','Select Input Images'));
elseif exist('ImgVols','var')
    save
    for i=1:length(ImgVols)
        V(i) = spm_vol(ImgVols(i).fname);
    end
    clear ImgVols;
    ImgVols=V;
    clear V;
end
if ~exist('optimal','var')
    optimal = spm_input('Normalize to Whole Head/GM/WM ?',1,'m','Normalize to Whole Head|Normalize to GM prior|Normalize to WM prior',[0 1 2], 1);
end
rotate = 0;%spm_input('Adjust pitch on oblique images (0.5 radians)?','+1','m','Yes|No',[1 0], 1);

if sge
    %     srccall =    'source /usr/local/sge/default/common/settings.sh';
    %     system(srccall);
    % Xvfb :9 -auth /etc/X99.cfg -screen scrn 800x600x8 &
    Xserver = spm_input('Xserver to display graphics to?','+1','s',['mr-bill:1.0']);
    stime=[date,'_',num2str(now)];
    results_dir = spm_input('Directory to save results in?','+1','s',[pwd,filesep,'template_',stime])
    cd('/')
    if ~exist(results_dir)
        mkdir(results_dir)
    end
    cd(results_dir)
    if ~exist('ImgVols','var')
        ImgVols = spm_vol(spm_get(Inf,'*.IMAGE','Select Input Images'));
    end
    %rotate = spm_input('Adjust pitch on oblique images (0.5 radians)?','+1','m','Yes|No',[1 0], 1);
    
    if rotate
      ImgVols=adjust_pitch(ImgVols,0.5);
    end
    sge_make_template(modality,ImgVols,results_dir,stime,optimal,Xserver);
else
  
  if ~exist('ImgVols','var')
    ImgVols = spm_vol(spm_get(Inf,'*.IMAGE','Select Input Images'));
  end
  if rotate
    ImgVols=adjust_pitch(ImgVols,0.5);
  end
  if strcmpi(modality,'t1')
    MNItemplate{1} = fullfile(spm('Dir'),'templates','T1.mnc')
  elseif strcmpi(modality,'t2')
    MNItemplate{1} = fullfile(spm('Dir'),'templates','T2.mnc');
  elseif strcmpi(modality,'pet')
    MNItemplate{1} = fullfile(spm('Dir'),'templates','PET.mnc');
  elseif strcmpi(modality,'pd')
    MNItemplate{1} = fullfile(spm('Dir'),'templates','PD.mnc');
  elseif strcmpi(modality,'epi')
    MNItemplate{1} = fullfile(spm('Dir'),'templates','EPI.mnc');
  elseif strcmpi(modality,'spect')
    MNItemplate{1} = fullfile(spm('Dir'),'templates','SPECT.mnc');
  elseif (strcmpi(modality,'transm') | strcmpi(modality,'other'))
    MNItemplate{1} = fullfile(spm('Dir'),'templates','Transm.mnc');
  else
    [template,priors]=make_template(parr,modality,ImgVols);
  end
  SegTemplateWeight=1; %get_weight(1,'SEG template');
  NrmTemplateWeight=1; %get_weight(1,'NRM template');
  write_seg=0;%input('save intermediate imgs? [0/1]');
  [template, priors, imgVols]=make_template_auto(ImgVols,MNItemplate,1,1,0,optimal,0,NrmTemplateWeight,0,SegTemplateWeight,0,write_seg);
  %[V,VG0,graphicmode,estimate,modulation,optimal,NRMsourceweight, NRMtemplateweight,SEGsourceweight, SEGtemplateweight,ApplyMask,write_seg] = deal(varargin{1:nargin});
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 sge_make_template(modality,ImgVols,results_dir,stime,optimal,Xserver)
make_template_auto_manager(ImgVols,modality,1,1,0,optimal,0,1,0,1,0,0,results_dir,stime,Xserver);
