Description of remove

0001 function remove_physio(proc,sinfo)
0002 % remove_physio.m
0003 %
0004 % Removes variance associated with cardiac and respiration regressors on a
0005 % slice-by-slice basis.
0006 %
0007 % Analysis steps to perform are specified in params which is a cell array of
0008 % structures.
0009 %
0010 % proc{1}.type  {'eval_model','calc_stats','spat_reg'}
0011 %
0012 % Different analysis steps require additional control parameters which are
0013 % specified in fields in the structure:
0014 %  eval_model:
0015 %    add_mean_to_resid - % Add the mean back to the estimated residuals.
0016 %    conlist
0017 %    Fconlist
0018 
0019 % 01/08/06 Petr Janata - apdated into function from remove_physio used for basic_states project
0020 
0021 % Handle input arguments.  We need at least analysis parameters and subject
0022 % info.
0023 error(nargchk(2,2,nargin,'struct'))
0024 
0025 % Read in the project globals
0026 % bs_globals
0027 
0028 nproc = length(proc);
0029 defidx = [];
0030 for iproc = 1:nproc
0031   if strcmp(proc{iproc}.type,'defs')
0032     defidx = iproc;
0033     rootdir = proc{iproc}.rootdir;
0034     verbose = proc{iproc}.verbose;
0035   end
0036 end
0037 
0038 stats_dir = 'stats';  % directory name within FEAT directory to store results in
0039 
0040 logstub = '';
0041 if ~isempty(logstub)
0042   logfile = fullfile(proc{defidx}.logdir, sprintf('%s_%s.txt',logstub,datestr(datenum(now),'yymmdd_HHMM')));
0043   logfid = fopen(logfile,'wt');
0044   fprintf('Writing logfile: %s\n', logfile);
0045 else
0046   logfid = 1;
0047 end
0048 
0049 if proc{defidx}.use_mc moco_str = '_mc'; else moco_str = ''; end
0050 
0051 start_dir = pwd;
0052 
0053 nsub_proc = length(sinfo);
0054 for isub = 1:nsub_proc
0055   subid = sinfo(isub).id;
0056   sdirs = check_subdirs(rootdir,subid,verbose);
0057 
0058   % Need to insert a session loop here
0059   
0060   % Load the physio information
0061   phys_fname = fullfile(sdirs.subdir, sprintf('%s_phys.mat', subid));
0062   fprintf(logfid,'Loading physio data from: %s\n', phys_fname);
0063   load(phys_fname)
0064   
0065   nruns = length(sinfo(isub).use_runs);
0066   for irun = 1:nruns
0067     run_id = sinfo(isub).use_runs(irun);
0068     
0069     % Make sure that run id in the phys structure matches the current run id
0070     if ~run_id == phys.ri(irun).id
0071       error(sprintf('Run ID mismatch: phys.ri.id=%d, run_id=%d', phys.ri(irun).id, run_id))
0072     end
0073     
0074     rdirs = check_rundirs(sdirs.physiodir, run_id,verbose);
0075     
0076     for iproc = 1:nproc
0077       cp = proc{iproc};
0078       switch cp.type
0079     case 'eval_model'
0080     
0081       epifname = fullfile(sdirs.epidir, sprintf('Y_run%d%s_bet.nii.gz', run_id, moco_str));
0082       % Check existence of the EPI file.  For some resting runs, data were not
0083       % originally converted, so we want to skip over those
0084       try
0085         check_exist(epifname);
0086       catch
0087         fprintf(logfid,'WARNING: Wanted run %d but could not find %s\n', run_id, epifname);
0088         continue
0089       end
0090     
0091       fsl_str = sprintf('avwval %s dim1', epifname);
0092       [status, nx] = unix(fsl_str);
0093       nx = str2num(nx);
0094 
0095       fsl_str = sprintf('avwval %s dim2', epifname);
0096       [status, ny] = unix(fsl_str);
0097       ny = str2num(ny);
0098 
0099       fsl_str = sprintf('avwval %s dim3', epifname);
0100       [status, nz] = unix(fsl_str);
0101       nslices = str2num(nz);
0102       
0103       % Estimate the threshold intensity that we want to use for this run and
0104       % this subject
0105       fsl_str = sprintf('avwstats %s -M', epifname);
0106       [status, mean_intensity] = unix(fsl_str);
0107       mean_intensity = str2num(mean_intensity);
0108       thresh = mean_intensity*cp.intensity_cutoff/100;
0109 
0110       % Calculate the mean of the EPI volume
0111       [fpath, fstub, fext] = fileparts(epifname);
0112       meanfname = fullfile(fpath, sprintf('%s_mean.nii.gz', strtok(fstub,'.')));
0113       fsl_str = sprintf('/usr/local/fsl/bin/avwmaths %s -Tmean %s', epifname, meanfname);
0114       fprintf(logfid,'%s\n', fsl_str);
0115       status = unix(fsl_str);
0116       if status
0117         error('Failed to calculate mean of input EPI data file')
0118       end
0119 
0120       for islice = 1:nslices
0121         % Specify the output directory
0122         featdir = fullfile(rdirs.rundir, sprintf('slice%02d.feat', islice));  % just the
0123         % stub.
0124         check_dir(featdir);
0125     
0126         funcfname = fullfile(sdirs.epidir, sprintf('slice%02d_data.nii.gz', islice));
0127     
0128         % Create the fsf file
0129         fsf = create_fsf;
0130         fsf = set_fsf_defaults(fsf, 'remove_physio');
0131     
0132         fsf.tr = cp.TR;
0133         fsf.fsldir = featdir;  % determines where design.fsf will be written
0134         fsf.outputdir = featdir;
0135         fsf.feat_files{1} = funcfname;
0136         fsf.npts = phys.ps.sinfo.nvol(run_id);
0137     
0138         % Add the physio regressors
0139         nsets = length(phys.ri(irun).ev_info);
0140         nev = 0;
0141         ev = create_ev;
0142         ev_names = {};
0143         
0144         ev_info = phys.ri(irun).ev_info;
0145         for iset = 1:nsets
0146           nev_in_set = ev_info(iset).nev;
0147           for iev = 1:nev_in_set
0148         cev = create_ev;
0149         cev.shape = 2;
0150         cev.fname = ev_info(iset).fnames{islice,iev};
0151         
0152         nev = nev+1;
0153         ev(nev) = cev;
0154         ev_names{nev} = ev_info(iset).ev_names{iev};
0155           end % for iev
0156         end % for iset
0157         fsf.ev = ev;
0158     
0159         fsf.evs_orig = nev;
0160         fsf.evs_real = nev;
0161 
0162         %
0163         % Deal with EV orthogonalization
0164         % At the very least, we have to initialize these to zero
0165         %
0166         for iev = 1:nev
0167           fsf.ev(iev).ortho = zeros(1,nev);
0168         end
0169     
0170         % Specify contrasts if we want to run the stats on the individual
0171         % regressors.
0172         %
0173         % NOTE: This has to be done at this stage rather than after packing all
0174         % of the slice data into a single volume so that the proper template files
0175         % are created for us to copy into the full volume analysis directory
0176 
0177         fsf = setup_fsl_con(fsf, ev_names, cp.conlist, cp.Fconlist);
0178 
0179         % Write the fsf file
0180         write_fsf(fsf);
0181     
0182         % Retain a copy of the fsf structure so that we can use it during
0183         % subsequent statistical evaluation
0184         fsf_fname = fullfile(featdir, 'fsf.mat');
0185         save(fsf_fname, 'fsf');
0186 
0187         % Switch to the directory that we'll be evaluating the model in
0188         cd(featdir)
0189     
0190         % Extract the relevant slice from the EPI data file
0191         fsl_str = sprintf('avwroi %s %s 0 %d 0 %d %d 1 0 %d', epifname, funcfname, nx, ny, islice-1, phys.ps.sinfo.nvol(run_id));
0192         fprintf(logfid,'%s\n', fsl_str);
0193         status = unix(fsl_str);
0194         if status
0195           error('Failed to extract relevant slice from the EPI data file')
0196         end
0197 
0198         % Extract the relevant slice from the mean data file
0199         slicemeanfname = fullfile(sdirs.epidir, sprintf('slice%02d_data_mean.nii.gz', islice));
0200         fsl_str = sprintf('avwroi %s %s 0 %d 0 %d %d 1 0 1', meanfname, slicemeanfname, nx, ny, islice-1);
0201         fprintf(logfid,'%s\n', fsl_str);
0202         status = unix(fsl_str);
0203         if status
0204           error('Failed to extract relevant slice from the mean data file')
0205         end
0206     
0207         % Evaluate the model.
0208         % Rather than calling FEAT which spawns a whole bunch of jobs, let's just
0209         % run those components that we really want
0210         
0211         % Check model integrity
0212         fsl_str = '/usr/local/fsl/bin/feat_model design';
0213         
0214         fprintf(logfid,'%s\n', fsl_str);
0215         status = unix(fsl_str);
0216         if status
0217           error('Checking of model failed')
0218         end
0219     
0220         % Remove the old stats directory
0221         if exist(stats_dir)
0222           unix_str = sprintf('rm -rf %s', stats_dir);
0223           unix(unix_str);
0224         end
0225         
0226         % Run the model
0227         fsl_str = sprintf(['/usr/local/fsl/bin/film_gls -rn %s -noest' ...
0228           ' %s design.mat %1.4f'], stats_dir, funcfname, thresh);
0229         fprintf(logfid,'%s\n', fsl_str);
0230         status = unix(fsl_str);
0231         if status
0232           error('FEAT run failed or was aborted')
0233         end
0234     
0235         if cp.add_mean_to_resid
0236           % Add the mean back into the residuals
0237           % NOTE: Adding the mean back has pros and cons.  It is good in the case that
0238           % the residuals are passed into routines that extract the brain based on a
0239           % proportion of the mean intensity in the image.  However, it is bad for some
0240           % routines (like my cardiac-related BOLD estimation) that use the residuals
0241           % image for masking.
0242 
0243           resid_fname = fullfile(stats_dir, 'res4d.nii.gz');
0244           fsl_str = sprintf('/usr/local/fsl/bin/avwmaths %s -add %s %s', resid_fname, slicemeanfname, resid_fname);
0245           fprintf(logfid,'%s\n', fsl_str);
0246           status = unix(fsl_str);
0247           if status
0248         error('Failed to add mean back into the residuals')
0249           end
0250         end % if cp.add_mean_to_resid
0251     
0252         % Remove the temporary functional data
0253         delete(funcfname);
0254         delete(slicemeanfname);
0255     
0256         cd(start_dir)
0257       end % for islice
0258 
0259       % Pack all of the slice by slice residuals and parameter estimates into
0260       % single volumes
0261             
0262       % Get a list of the files that we want to pack together
0263       stats_flist = dir(fullfile(rdirs.rundir,'slice01.feat/stats/*.nii.gz'));
0264       nstats_images = length(stats_flist);
0265       for iimg = 1:nstats_images
0266         merge_list = {};
0267         for islice = 1:nslices
0268           merge_list{islice} = fullfile(rdirs.rundir, sprintf('slice%02d.feat', islice),'stats',stats_flist(iimg).name);
0269         end
0270     
0271         outfname = fullfile(rdirs.fullvolstatsdir,stats_flist(iimg).name);
0272         fsl_str = sprintf('/usr/local/fsl/bin/avwmerge -z %s %s', outfname, cell2str(merge_list, ' '));
0273         fprintf(logfid,'%s\n', fsl_str);
0274         status = unix(fsl_str);
0275         if status
0276           error(sprintf('Failed to merge slice files into: %s', outfname))
0277         end
0278       end % for iimg = 1:nstats_images
0279       % END OF EVAL_MODEL
0280     
0281     case 'calc_stats'  % Calculate significance maps using FEAT
0282       fprintf(logfid, sprintf(['\nEvaluating statistical significance of the' ...
0283         ' physiological variables in run %d\n'], run_id));
0284       slice_featdir = fullfile(rdirs.rundir,'slice01.feat');
0285 
0286       % Copy the mask that is associated with the original EPI volume
0287       orig_maskfname = fullfile(sdirs.epidir, sprintf('Y_run%d%s_bet_mask.nii.gz', run_id, moco_str));
0288       maskfname = fullfile(fullvoldir,'mask.nii.gz');
0289       unix_str = sprintf('cp %s %s', orig_maskfname, maskfname);
0290       status = unix(unix_str);
0291       if status
0292         error('Error copying mask file')
0293       end
0294       
0295       % Copy the degrees of freedom file from first slice analysis
0296       unix_str = sprintf('cp %s %s', fullfile(slice_featdir,stats_dir,'dof'), rdirs.fullvolstatsdir);
0297       status = unix(unix_str);
0298       if status
0299         error('Error copying dof file')
0300       end
0301 
0302       % Copy the original full volume time series
0303       epifname = fullfile(sdirs.epidir, sprintf('Y_run%d%s_bet.nii.gz', run_id, moco_str));
0304       funcfname = fullfile(rdirs.fullvoldir,'filtered_func_data.nii.gz');
0305       unix_str = sprintf('cp %s %s', epifname, funcfname);
0306       fprintf(logfid,'%s\n', unix_str);
0307       status = unix(unix_str);
0308       if status
0309         error('Error copying original EPI data')
0310       end
0311     
0312       % Extract a single volume that we will use for overlays
0313       example_funcfname = fullfile(rdirs.fullvoldir, 'example_func.nii.gz');
0314       fsl_str = sprintf('/usr/local/fsl/bin/avwroi %s %s 0 1', funcfname, example_funcfname);
0315       fprintf(logfid,'%s\n', fsl_str);
0316       status = unix(fsl_str);
0317       if status
0318         error('Error extracting example functional volume')
0319       end
0320     
0321       % Load the fsf structure from slice 1 as a template
0322       fsf_fname = fullfile(slice_featdir,'fsf.mat');
0323       check_exist(fsf_fname);
0324       load(fsf_fname)
0325 
0326       % Modify the fsf structure accordingly
0327       fsf = set_fsf_defaults(fsf, 'evaluate_physio');
0328       fsf.outputdir = rdirs.fullvoldir;
0329       fsf.fsldir = rdirs.fullvoldir;
0330       fsf.feat_files{1} = rdirs.fullvoldir;  % FEAT directory
0331     
0332       write_fsf(fsf);
0333     
0334       % Save an fsf.mat file for posterity's sake
0335       save(fullfile(rdirs.fullvoldir, 'fsf.mat'), 'fsf')
0336 
0337       % Move to the 'FEAT' directory
0338       cd(rdirs.fullvoldir)
0339     
0340       fsl_str = sprintf('/usr/local/fsl/bin/feat %s', 'design.fsf');
0341       fprintf(logfid,'%s\n', fsl_str);
0342       status = unix(fsl_str);
0343       if status
0344         error('Problems running FEAT')
0345       end
0346       cd(start_dir)
0347       
0348       % end of calc_stats section
0349       
0350     case 'spat_reg'  % perform spatial registration
0351       % Load the existing fsf structure
0352       fsf_fname = fullfile(rdirs.fullvoldir,'fsf.mat');
0353       check_exist(fsf_fname);
0354       load(fsf_fname)
0355 
0356       fsf = set_fsf_defaults(fsf, 'coreg');
0357       fsf.outputdir = rdirs.fullvoldir;
0358       fsf.fsldir = rdirs.fullvoldir;
0359       if ~iscell(fsf.feat_files)
0360         fsf.feat_files = {rdirs.fullvoldir};
0361       end
0362       
0363       % Specify coplanar
0364       coplanar_stub = sprintf('%s_coplanar', subid);
0365       coplanar_fname = fullfile(sdirs.anatdir,sprintf('%s.nii.gz', coplanar_stub));
0366       fsf.initial_highres_files{1} = coplanar_fname;
0367       
0368       % Specify hires
0369       hires_stub = sprintf('%s_hires', subid);
0370       hires_fname = fullfile(sdirs.anatdir,sprintf('%s.nii.gz', hires_stub));
0371       fsf.highres_files{1} = hires_fname;
0372       
0373       write_fsf(fsf);
0374 
0375       % Move to the 'FEAT' directory
0376       cd(rdirs.fullvoldir)
0377       
0378       fsl_str = sprintf('/usr/local/fsl/bin/feat %s', 'design.fsf');
0379       fprintf(logfid,'%s\n', fsl_str);
0380       status = unix(fsl_str);
0381       if status
0382         error('Problems running FEAT')
0383       end
0384       
0385       % Save an fsf.mat file for posterity's sake
0386       save(fullfile(rdirs.fullvoldir, 'fsf.mat'), 'fsf')
0387 
0388       cd(start_dir)
0389       % end of spat_reg block
0390     otherwise
0391       end % switch proc.type
0392     end % for iproc=
0393   end % for irun=
0394 end % for isub=
0395 end % remove_physio
0396 
0397 function sdirs = check_subdirs(rootdir,subid, verbose)
0398   subdir = fullfile(rootdir, subid);
0399   epidir = fullfile(subdir, 'epi');
0400 
0401   fsldir = fullfile(subdir, 'fsl');
0402   check_dir(fsldir, verbose);
0403 
0404   anatdir = fullfile(subdir, 'anat');
0405   check_dir(anatdir);
0406   
0407   physiodir = fullfile(subdir, 'physio_resid');
0408   check_dir(physiodir, verbose);
0409 
0410   % Populate the output struct
0411   sdirs = struct('subdir', subdir, 'epidir', epidir, 'fsldir', fsldir, ...
0412       'anatdir', anatdir, 'physiodir', physiodir);
0413 end % check_subdirs
0414 
0415 function rdirs = check_rundirs(rootdir, runid, verbose)
0416   rundir = fullfile(rootdir,sprintf('run%d', runid));
0417   check_dir(rundir, verbose);
0418   
0419   fullvoldir = fullfile(rundir, 'fullvol');
0420   check_dir(fullvoldir, verbose);
0421     
0422   fullvolstatsdir = fullfile(fullvoldir, 'stats');
0423   check_dir(fullvolstatsdir, verbose);
0424   
0425   % Populate the output struct
0426   rdirs = struct('rundir', rundir, 'fullvoldir', fullvoldir, ...
0427       'fullvolstatsdir', fullvolstatsdir);
0428 end % check_rundirs
remove_physio

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
