function vbm_custom3(runmode,templatedir,optimal,Anatnames,other_img_dir,results_dir,sm)
%VBM_custom3 plus other images ...
% this routine performs VBM preprocessing steps on a number of T1 images.
% The user is prompted for the images, as well as custom template and
% priors.  The normalization is performed (for now) either GM -> GM, or
% WM -> WM.  i.e., each subjects WM or GM is normalised to the WM or GM
% apriori image.
% Additionally, the resulting DCT deformation is applied to related images
% for each subject.  These images could be, e.g., Leukariosis masks from
% HistoSeg output, or they could be Log(Jacobian) maps from a serial fluid
% warping for each subject.  These masks are normalized into [1 1 1] voxel
% size, with NN interpolation.... those settings could be changed if you
% look for them in the code.
% USAGE: 
%   vbm_custom3  prompts user for all inputs.
%   vbm_custom3('sge')  runs over grid engine.
%   vbm_custom3(runmode,templatedir,optimal,AnatNames,other_img_dir,results_dir,sm)
%   INPUTS:
%   runmode = 'sge' run over the Grid.
%             'uni' run on this computer only.
%   templatedir = '/my/template/directory' - this directory must contain
%   s8template, s8gray, s8white, and s8csf (.img, .hdr pairs)
%   optimal = 0 : normalize whole head to whole head template
%             1 : normalize GM segmentation to GM prior
%             2 : normalize WM segmentation to WM prior
%
%   AnatNames = cell array of absolute filenames of input images
%   other_img_dir = '/my/other/dir' - this directory should contain extra
%   images to be warped for each subject.
%   results_dir = '/my/results/dir'  - directory to store results in.
%   sm = FWHM(mm) smoothing to apply to GM images [0=none]
%
% SEE ALSO: vbm_custom, vbm_custom2
% Matt Senjem
spm_defaults;
global defaults
if ~exist('runmode','var')
    runmode = 'uni';
end
switch runmode
    case 'sge'
        sge=1;
    otherwise
        sge=0;
end
if exist('templatedir','var')
    template_cell{1} = [templatedir,'/s8template.img'];
    template_cell{2} = [templatedir,'/s8gray.img'];
    template_cell{3} = [templatedir,'/s8white.img']; 
    template_cell{4} = [templatedir,'/s8csf.img'];
end
%optimal = spm_input('Normalize to GM/WM prior?',1,'m','Normalize to GM prior|Normalize to WM prior',[1 2], 1);
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
if ~exist('Anatnames','var')
    Anatnames = spm_get(Inf,'*.IMAGE',{'Select Anat Images'});
end
Anatdir  = get_rootdir(Anatnames);
slash = findstr(filesep,Anatdir);
rootdir = Anatdir;%(1:slash(end)-1);

if ~exist('other_img_dir','var')
    other_img_dir  = spm_input('Other img rootdir?','+1','s',rootdir);
end

if ~exist('results_dir','var')
    results_dir = spm_input('Directory to save results in?','+1','s',rootdir)
end

if ~exist('email','var')
    email = spm_input('Email address to notify?','+1','s','Senjem.Matthew1@mayo.edu');
end
if ~exist('priority','var')
    priority = spm_input('Job priority?','+1','i',1);
end
priority = max(-1024,priority);
priority = min(1024,priority);

if ~exist('sm','var')
    sm = spm_input('FWHM(mm) smoothing [0=none]?','+1','w','12');
end
if ~exist('stime','var')
    stime=[date,'_',num2str(now)];%stime = date;
end
other_img_names = find_AnatJacpairs(Anatnames,other_img_dir);
cd('/')
if ~exist(results_dir)
    mkdir(results_dir)
end

disp(['pwd = ',pwd])
cd(results_dir)
disp(['pwd = ',pwd])

numsubjects = length(Anatnames);

if ~exist('template_cell','var')
    template_cell{1} = spm_get(1,'s8template*.IMAGE','Select template for normalization'); 
    template_cell{2} = spm_get(1,'s8gray*.IMAGE','Select gray matter prior'); 
    template_cell{3} = spm_get(1,'s8white*.IMAGE','Select white matter prior'); 
    template_cell{4} = spm_get(1,'s8csf*.IMAGE','Select CSF prior'); 
end


if sge
    save vbm_custom3.mat
    sge_mls_fluid_prep2('norm2custom',template_cell,'smooth',sm,...
        'Anatnames',Anatnames,'logJnames',other_img_names,...
        'results_dir',results_dir,'stime',stime,'optimal',optimal,'email',email,'priority',priority);
else
save
disp('not yet implemented in single-CPU mode...')
disp('try vbm_custom3(''sge'')')
end

return

function other_img_names = find_AnatJacpairs(Anatnames,other_img_dir,other_img_names)
subdirs = dir(other_img_dir);
if ~exist('other_img_names','var')
    other_img_names = {};
end
for k=1:length(Anatnames)
    fnd=0;
    
    disp(['matching ',Anatnames{k}])
    [pth,nam,ext] = fileparts(Anatnames{k});
    uscore = findstr('_',pth);
    if ~isempty(uscore)
        namparts = get_name_parts(nam);
        disp(char(namparts))
        for n=1:length(namparts)
            while strcmp(namparts{n}(1),'0')
                namparts{n}=namparts{n}(2:end);
            end
            if fnd<2
                disp(' ')
                disp(' ')
                disp(['searching : ',fullfile(other_img_dir,['*',namparts{n},'*',ext])])
                files = dir(fullfile(other_img_dir,['*',namparts{n},'*',ext]));
                disp(' ')
                i22=0;
                for i11=1:length(files)
                    disp(files(i11).name)
                    if ~strcmp(files(i11).name,[nam,ext])
                        files2(i11-i22)=files(i11);;
                    else
                        i22=i22+1;
                    end
                end
                clear files
                files=files2;
                clear files2;
                disp(' ')
                disp('found some matching files ...')
                disp('---------------------------------------------------------------')
                disp(char(files.name))
                keep = 'a';%input('warp these img files? ([a]ll /[n]one/[s]ome)','s');
                if isempty(keep) | strcmpi(keep(1),'a')
                    fnd=2;
                    for ppp=1:length(files)
                        [other_img_names{k,ppp}] = fullfile(other_img_dir,files(ppp).name);
                    end
                    fnd=2;
                    disp(' ')
                    disp('---------------------------------------------------------------')
                    disp(char(other_img_names{k,:}))
                    disp(' ')
                    break
                elseif strcmpi(keep(1),'n')
                    continue
                elseif strcmpi(keep(1),'s')
                    disp('searching ... ')
                else
                    error('incorrect input ... sorry I don''t know know what to do next')
                end
                disp(' ')
                disp('---------------------------------------------------------------')
            end
            if fnd<2
                for i1=1:length(files)
                    if ~strcmp(files(i1).name,[nam,ext])
                        disp(' ')
                        disp('---------------------------------------------------------------')
                        disp(['found: ',files(i1).name])
                        disp(' ')
                        keep = input('warp this img? ([y]es/[n]o/[d]one)','s');
                        if strcmpi(keep(1),'y')
                            fnd=1;
                            other_img_names{k,end+1} = fullfile(other_img_dir,files(i1).name);
                        elseif strcmpi(keep(1),'d')
                            fnd=2;
                            disp(' ')
                            disp('---------------------------------------------------------------')
                            disp(char(other_img_names{k,:}))
                            disp(' ')
                            
                            
                            break
                        elseif strcmpi(keep(1),'n')
                            disp('skipping ... ')
                        else
                            error('incorrect input ... sorry I don''t know know what to do next')
                        end
                    end
                end
            end
        end
    end    
    if fnd<2
        for i2=1:length(subdirs)
            if subdirs(i2).isdir  & ~strcmp('.',subdirs(i2).name) & ~strcmp('..',subdirs(i2).name)
                disp(['scanning ',fullfile(other_img_dir,subdirs(i2).name),' ... '])
                other_img_names{k,end+1} = find_AnatJacpairs({Anatnames{k}},[other_img_dir,filesep,subdirs(i2).name],other_img_names);
            end
        end
    end
end
return

function parts = get_name_parts(nam)
uscore = findstr('_',nam);
parts{1}=nam(1:uscore(1)-1);
if length(uscore)>1
    for i2=1:length(uscore)-1
        parts{i2+1} = nam(uscore(i2)+1:uscore(i2+1)-1);
    end
    parts{i2+2}=nam(uscore(i2+1)+1:end);
end
order=1:length(parts);
maxlen=length(parts{1});
% for p=1:length(parts)
%     plen(p)=length(parts{p});
% end
% [a,b]=sort(plen);
% [parts{1:length(parts)}] = deal(parts{fliplr(b)});
return
