spike_check

  status = spike_check(sinfo, inpath, outpath, scan_offset, nvox, spkcrit);

  Checks GE format data for spikes and dropouts in the signal.

  nvox -- size of volume in voxels (x,y,z) centered on original volume to use
  in the analysis
 
  spkcrit -- a structure containing fields with criteria for calling something
  a spike:
     .nstd - if the difference between two samples is this many standard
             deviations (computed across the single run), the sample will be
             flagged as bad. Default = 1.


  The script is based on Souheil Inati's QA.m script, but modified to handle
  multiple runs, and place various data in the output structure

0001 function [info, status] = spike_check(sinfo, inpath, outpath, scan_offset, nvox, sc) 
0002 %  status = spike_check(sinfo, inpath, outpath, scan_offset, nvox, spkcrit);
0003 %
0004 %  Checks GE format data for spikes and dropouts in the signal.
0005 %
0006 %  nvox -- size of volume in voxels (x,y,z) centered on original volume to use
0007 %  in the analysis
0008 %
0009 %  spkcrit -- a structure containing fields with criteria for calling something
0010 %  a spike:
0011 %     .nstd - if the difference between two samples is this many standard
0012 %             deviations (computed across the single run), the sample will be
0013 %             flagged as bad. Default = 1.
0014 %
0015 %
0016 %  The script is based on Souheil Inati's QA.m script, but modified to handle
0017 %  multiple runs, and place various data in the output structure
0018 %
0019 
0020 % 04/07/01 Petr Janata
0021 % 04/27/01 PJ Added ability to handle EPI runs distributed across multiple
0022 %             series in a single exam. To take advantage of this specification,
0023 %             the epi run numbers associated with a series must be specified in
0024 %             a 3rd column of cells in the series_mapping field of the subject
0025 %             information structure.  Otherwise, the cond_order field is used
0026 %             and assumes that all EPI runs in an exam were collected under a
0027 %             single series.
0028 % 04/29/01 PJ fixed indexing problem when grabbing number of volumes for a
0029 %             run. Needed to refer to cond_order field.
0030 % 01/30/01 PJ Brought in line with convert_format.m (no need to check cond_order)
0031 
0032 status = 0;
0033 
0034 nsub = length(sinfo);
0035 
0036 if nargin < 2
0037   inpath = './';
0038   outpath = './';
0039 elseif nargin < 3
0040   outpath = './';
0041 end
0042 
0043 if nargin < 4
0044   scan_offset = 0;
0045 end
0046 
0047 if nargin < 5
0048   nvox = [15 15 2];
0049 end
0050 
0051 if isempty(outpath)
0052   outpath = './';
0053 end
0054 
0055 if isempty(inpath)
0056   inpath = './';
0057 end
0058 
0059 if nargin < 6
0060   sc.nstd = 1;
0061 end
0062 
0063 for isub = 1:nsub
0064   
0065   s = sinfo(isub);
0066   subj_root = s.id;
0067   
0068   info(isub).id = s.id;
0069   
0070   num_epi = 0;
0071   for iexam = 1:s.nexams
0072 
0073     exam_root = sprintf('%05d', s.exam_nums(iexam));
0074     datapath = fullfile(inpath, subj_root, exam_root);
0075     
0076     series=s.series_mappings{iexam};
0077     nseries= size(series,1);
0078 
0079     for mapping_idx = 1:nseries
0080       indir = fullfile(datapath, char(series(mapping_idx,1)));
0081       
0082       switch char(series(mapping_idx,2))
0083     case {'epi','epi1','epi2','epi3','epi4','epi_12','epi_34'}
0084 
0085       % Check to see if there is condition order information in the series
0086       % mappings
0087       if size(series,2) == 3
0088         nruns = length(series{mapping_idx,3});
0089         multi_epi_series = 1;
0090       else
0091         nruns = length(s.cond_order);
0092         multi_epi_series = 0;
0093       end
0094 
0095       % Create the name of the first file in inDir
0096       firstfile = fullfile(indir,'I.001');
0097 
0098       % Read the header
0099       [su_hdr,ex_hdr,se_hdr,im_hdr,pix_hdr,im_offset] = GE_readHeader(firstfile);
0100 
0101       nX =  im_hdr.imatrix_X; % X Voxels
0102       nY =  im_hdr.imatrix_Y; % Y Voxels
0103       nZ =  im_hdr.slquant;   % Z Voxels
0104       volSize = [nX nY nZ];
0105 
0106       % Initialize
0107       imageVol = zeros(nX,nY,nZ);
0108       
0109       xa = floor(nX/2-nvox(1));
0110       xb = ceil(nX/2+nvox(1));
0111       ya = floor(nY/2-nvox(2));
0112       yb = ceil(nY/2+nvox(2));
0113       za = floor(nZ/2-nvox(3));
0114       zb = ceil(nZ/2+nvox(3));
0115       nx = xb-xa+1;
0116       ny = yb-ya+1;
0117       nz = zb-za+1;
0118       subVol = zeros(nx,ny,nz);
0119 
0120       % Check for number of EPI runs
0121       for irun = 1:nruns
0122         num_epi = num_epi+1;
0123         
0124         if multi_epi_series
0125           nvol = s.nvol(series{mapping_idx,3}(irun));
0126         else
0127           nvol = s.nvol(irun);
0128         end
0129         
0130         info(isub).exam{num_epi} = exam_root;
0131         info(isub).series{num_epi} = char(series(mapping_idx,1));
0132         info(isub).run(num_epi) = num_epi;
0133         
0134         info(isub).data{num_epi} = zeros(nvol,1);
0135         
0136         if multi_epi_series
0137           run_offset = sum(s.nvol(series{mapping_idx,3}(1:irun))) - s.nvol(series{mapping_idx,3}(irun));
0138         else
0139           run_offset = sum(s.nvol(1:irun))-s.nvol(irun);
0140         end
0141         
0142         start_num = run_offset+1+scan_offset; % increment 1 to get to first image
0143 
0144         for ivol = 1:nvol
0145           [imageVol, lastfile] = GE_readVolume(indir, start_num + ivol, volSize, 16, ...
0146           im_offset);
0147           subVol = imageVol(xa:xb,ya:yb,za:zb);
0148           info(isub).data{num_epi}(ivol) = mean(subVol(:));
0149           fprintf('runnum = %d\tvolnum = %d\tread %s\n',num_epi, ivol,lastfile);
0150         
0151         end % for ivol=
0152         
0153         %
0154         % Calculate some stats
0155         %
0156         tmp = info(isub).data{num_epi};
0157         info(isub).mean(num_epi) = mean(tmp);
0158         info(isub).std(num_epi) = std(tmp);
0159         
0160         % Find locations where there are sharp transients
0161         tmpdiff = abs(diff(info(isub).data{num_epi})); 
0162         badidx = find(tmpdiff > sc.nstd * info(isub).std(num_epi));
0163         info(isub).badidx{num_epi} = badidx;
0164         
0165         info(isub).numbad(num_epi) = fix(length(badidx)/2); % **rough**
0166                                                                 % estimate of
0167                                                                 % number of
0168                                                                 % spikes
0169                                 
0170       
0171       end % for irun=
0172       end % switch
0173     end % for mapping_idx =
0174   end % for iexam =
0175 end % for isub=