function varargout=sge_mls_fluid_prep2(varargin)
% This is a big, nasty routine to do fMRI SPM analysis
% Generally,
% To analyze BRIK files from a single subject, you should use fmri_new_spm.
% Not to worry, you WILL be prompted for any missing input arguments.
%
% this function will:
%      1. convert AFNI (*.BRIK, *.HEAD) files to AnalyzeImage(7.5) (*.img , *.hdr)
%         using the perl script brik2img1 installed on mr-matrix.mayo.edu in /usr/local/bin/
%      2. Try to distinguish Anatomical SPGR images from functional EPI images. (usually works)
%      3. Flip and re-orient the images in such that they are appropriate for SPM analysis.
%      4. co-register EPI volumes of each subject to the first EPI for that subject, creating a mean*.img
%      5. Use NMI to register mean*.img to the Anatomical scan , and apply these parameters to all EPI images.
%      6. (OPTIONAL) create a custom template from all the SPGR images.
%      7. Normalize T1 images to custom or MNI template, & apply these parameters to EPI images.
%      8. Smooth the Normalized EPI volumes by convolution with a Gaussian filter of the specified FWHM (defaults to 9 mm)
%      9. Enter the smoothed volumes into an SPM.mat data structure, and analyze the model,
%         complete with individualized contrasts for each session/subject.
% To create a custom template and normalize to it, using BRIK input
% files:
% mls_fmri_prep('BRIK','numsubjects',5,'sessionspersubject',3,'scanspersession',64,'smooth',9,'mktemplate','norm2custom')
% mls_fmri_prep('fmri','numsubjects',1,'sessionspersubject',5,'scanspersession',64,'stats')
% To create a custom template and normalize to it:
% mls_fmri_prep('numsubjects',5,'sessionspersubject',3,'scanspersession',64,'smooth',9,'mktemplate','norm2custom')
%
% To prompt the user for a template & priors to which to normalize:
% mls_fmri_prep('numsubjects',5,'sessionspersubject',3,'scanspersession',64,'smooth',9,'norm2custom')
%
% To normalize to MNI template
% mls_fmri_prep('numsubjects',5,'sessionspersubject',3,'scanspersession',64,'smooth',9,'norm2mni')
%
% Author: Matt Senjem, May 05, 2004.
% ver. 1.0
% Written by : Matt Senjem May 5, 2004.
%______________________________________________________________________________________________%
% ______________________________________________________________________________________________%

spm_defaults;
global defaults;
if nargin<1
    disp('Warning: no inputs specified !')
end
vbmLevel=0;
global defaults
if nargin==1 & iscell(varargin)
    foo=varargin{1};
    clear nargin
    narg=length(foo);
    clear varargin;
    varargin=foo;
    clear foo;
else
    narg=nargin;
end
k=1;
while k<=narg
    if isstr(varargin{k}) & strcmpi(varargin{k},'numsubjects')
        if isnumeric(varargin{k+1})
            k=k+1;
            numsubjects=varargin{k};
        end
    elseif isstr(varargin{k}) & strcmpi(varargin{k},'stime')
        if ischar(varargin{k+1})
            k=k+1;
            stime = varargin{k};
        end
    elseif isstr(varargin{k}) & strcmpi(varargin{k},'optimal')
      if isnumeric(varargin{k+1})
          k=k+1;
          optimal = varargin{k};
      end
    elseif isstr(varargin{k}) & strcmp(varargin{k},'Anatnames')
        if iscell(varargin{k+1})
            k=k+1;
            Anatnames=varargin{k};
        end
    elseif isstr(varargin{k}) & strcmpi(varargin{k},'email')
        if isstr(varargin{k+1})
            k=k+1;
            email=varargin{k};
        end
    elseif isstr(varargin{k}) & strcmpi(varargin{k},'priority')
        if isnumeric(varargin{k+1})
            k=k+1;
            priority=varargin{k};
        end
    
    elseif isstr(varargin{k}) & strcmp(varargin{k},'logJnames')
        if iscell(varargin{k+1})
            k=k+1;
            logJnames=varargin{k};
        end
    elseif isstr(varargin{k}) & strcmpi(varargin{k},'fluid')
        modality='fluid';
    elseif isstr(varargin{k}) & strcmpi(varargin{k},'norm2custom')
        norm2custom=1;
        a=varargin{k+1};
        if iscell(varargin{k+1})
            ctemplate_cell =varargin{k+1};
            k=k+1;
        end
    elseif isstr(varargin{k}) & strcmpi(varargin{k},'norm2mni')
        norm2mni=1;
        if isstr(varargin{k+1})
            if strcmpi(varargin{k+1},'t1')
                MNItemplate  = spm_vol(fullfile(spm('Dir'),'templates','T1.mnc'))
                k=k+1;
            elseif strcmpi(varargin{k},'t2')
                MNItemplate  = spm_vol(fullfile(spm('Dir'),'templates','T2.mnc'));
                k=k+1;
            elseif strcmpi(varargin{k},'pd')
                MNItemplate  = spm_vol(fullfile(spm('Dir'),'templates','PD.mnc'));
                k=k+1;
            elseif strcmpi(varargin{k},'spect')
                MNItemplate  = spm_vol(fullfile(spm('Dir'),'templates','SPECT.mnc'));
                k=k+1;
            elseif strcmpi(varargin{k},'pet')
                MNItemplate  = spm_vol(fullfile(spm('Dir'),'templates','PET.mnc'));
                k=k+1;
            elseif strcmpi(varargin{k},'epi')
                MNItemplate  = spm_vol(fullfile(spm('Dir'),'templates','EPI.mnc'));
                k=k+1;    
            elseif strcmpi(varargin{k},'transm')
                MNItemplate  = spm_vol(fullfile(spm('Dir'),'templates','transm.mnc'));
                k=k+1;    
            else
                error('unknown modality')
            end
            ctemplate=MNItemplate;
        end
    elseif isstr(varargin{k}) & strcmpi(varargin{k},'smooth')
        smooth=1;
        if isnumeric(varargin{k+1})
            k=k+1;
            if varargin{k}~=0
                fwhm = [varargin{k} varargin{k} varargin{k} ];
            else
                smooth=0;
                fwhm=0;
            end
        end
    elseif isstr(varargin{k}) & strcmpi(varargin{k},'prefix')
        if isstr(varargin{k+1})
            k=k+1;
            prefx=varargin{k};
        end
    elseif isstr(varargin{k}) & strcmpi(varargin{k},'results_dir')
        if isstr(varargin{k+1})
            k=k+1;
            results_dir=varargin{k};
        end
    else
        save
        error(['unknown input ', varargin{k}])
    end
    k=k+1;
end
if ~exist('numsubjects','var')
  numsubjects = length(Anatnames);
end
if ~exist('stime','var')
    stime=[date,'_',num2str(now)];%stime = date;
end
modality='fluid';
cd('/')

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 = min(1024,max(-1024,priority));

if ~exist('results_dir','var')
    results_dir = spm_input('Directory to save results in?','+1','s',[pwd,filesep,stime]);
    if ~exist(results_dir,'dir')
        mkdir(results_dir);
    end
    
end
cd(results_dir)
disp(['pwd = ',pwd])
if (exist('ctemplate_cell','var') & length(ctemplate_cell)==4)
    ctemplate = spm_vol(ctemplate_cell{1})
    cpriors = char(ctemplate_cell{2:4})
end
if (exist('mktemplate','var') | exist('norm2mni','var'))
    if ~exist('MNItemplate','var')
        MNItemplate  = spm_vol(spm_get(Inf,'*.IMAGE','Select MNI template for normalization'));
    end
end
if exist('norm2custom','var') | exist('norm2mni','var')
    if exist('norm2custom','var') & norm2custom
        if ~exist('ctemplate','var')
            ctemplate = spm_vol(spm_get(1,'*.IMAGE','Select custom ctemplate for normalization'));
        end
        if ~exist('cpriors','var')
            priors = spm_get(3,'s*.IMAGE','Select gray/white and CSF prior images');
        else
            priors=cpriors;
        end

        if ~exist('NRMtemplateweight','var')
            NRMtemplateweight=0;
        end
        if ~exist('SEGtemplateweight','var')
            SEGtemplateweight=0;
        end
        if ~exist('bmask','var')
            bmask='';%fullfile(spm('Dir'),'apriori','brainmask.mnc');	%
        end
    elseif exist('norm2mni','var') & norm2mni
        if ~exist('dseg','var'),dseg=defaults.segment;end
        priors=      dseg.estimate.priors;
        NRMtemplateweight=1;
        SEGtemplateweight=1;
        bmask=fullfile(spm('Dir'),'apriori','brainmask.mnc');	% always use brainmask with MNI template

    end
end

if exist('Anatnames','var')
    rotate = 0;%spm_input('Adjust pitch on oblique images (0.5 radians)?','+1','m','Yes|No',[1 0], 1);
    imlist = Anatnames;
    for v=1:length(Anatnames)
        imgtmp =dir(imlist{v});
        imglist(v)=imgtmp;
    end
    AnatVols = spm_vol(char(Anatnames));
    if rotate==1    
        AnatVols=adjust_pitch(AnatVols,.5)
    end
    clear Anatnames;
end
%try
%    save
%    for i=1:length(AnatVols)
%        for j=1:size(char(logJnames{i,:}),1)
%            if ~isempty(logJnames{i,j})
%                EPIVols{i,j}   = spm_vol(logJnames{i,j});
%%            end
%        end
%    end
%catch 
%    try

savefile = ['MFPparams',stime,'.mat'];
save(savefile)
disp(['saved parameters to ',savefile]);
[s1,s2] = size(logJnames);
lan = length(AnatVols);
if s1==lan
    for i=1:lan
        %save
        for j=1:size(char(logJnames{i,:}),1)
            if ~isempty(logJnames{i,j})
                EPIVols{i,j}   = spm_vol(logJnames{i,j});
            end
        end
    end
    
elseif s2 == lan
    for i=1:lan
        %save
        for j=1:size(char(logJnames{:,i}),1)
            if ~isempty(logJnames{j,i})
                EPIVols{j,i}   = spm_vol(logJnames{j,i});
            end
        end
    end
    
else
    error('foo')
end
%catch
%   save
%       error('try loading matlab.mat')
%   end
%end

save(savefile)
disp(['saved parameters to ',savefile]);
% ______________________________________________________________________________________________%

% create a custom template
if exist('mktemplate','var') & mktemplate &  exist('AnatVols','var')
    [ctemplate ,cpriors] = make_template_auto(AnatVols,MNItemplate,1,1,0,1,0,1,0,1,0,0)
end

% ______________________________________________________________________________________________%
% normalize Anatomical to template (MNI or custom)
% and optionally apply the same parameters to the EPI / pet / Jacobian images.
if exist('norm2custom','var') |  exist('norm2mni','var')
    if ~exist('dseg','var'),dseg=defaults.segment;end

    for i=1:length(AnatVols)
        [pth,nam,ext] = fileparts(AnatVols(i).fname);
        matname0 = fullfile(pth,[nam '_sn.mat']);
        if exist(matname0,'file')
            matfiles{i}=matname0;
        end
    end
    if ~exist('matfiles','var')
        matfiles   = {};%spm_get([0 length(AnatVols)],'*_sn.mat',{'Select *_sn.mat files'});
    end
    %if ~exist('EPIVols','var')
    %    disp([])
    %    EPIVols=[];
    %end
    % dseg.estimate.priors=priors;
    % we should have the template & priors selected by now (I hope!!!)
    %% FILL HERE: DONE!! (beta)
    %if strcmpi(modality,'pet')
    %    voxsize=1;
    %elseif strcmpi(modality,'fluid')
        voxsize=1;
        %end
%try,    
  norm_fluid_auto_manager(AnatVols,matfiles,EPIVols,ctemplate,priors,1,1,...
        optimal,0,NRMtemplateweight,0,...
        SEGtemplateweight,0,2,1, ...
        0,1,bmask,voxsize,smooth,fwhm,results_dir,stime,email,priority);
%catch, save mfpERROR.mat,end

    %    norm_fmri_auto(AnatVols,matfiles,EPIVols,template,priors,graphicmode,modulation,...
    %		   optimal,NRMsourceweight,NRMtemplateweight,SEGsourceweight,...
    %		   SEGtemplateweight,ApplyMask,max_seg,writeseg,writeanat,writeepi,bmask);
    epj=0;
    for ep=1:length(AnatVols)
        for ej=1:size(EPIVols,2)
            try 
                if ~isempty(EPIVols{ep,ej})
                    epj=epj+1;
                    [Epth,Enam,Eext]=fileparts(EPIVols{ep,ej}.fname);
                    foo=fullfile(Epth,['w_' Enam Eext]);
                    epiNames{epj}=foo;
                    clear foo;
                    disp(['look for: ',epiNames{epj}])
                end
            catch
                save(['err_',mfilename,'_',num2str(epj)])
            end
        end
    end
end

