fmriquality

As said on the quality assessment page, this function is used to determine whether or not a given dataset (run) is suitable for subsequent data analysis.

To run the fMRI quality checking function, the images need to be in one of the functional imaging data formats currently supported by the xff class:

• BrainVoyager's FMR (STC) format
• BrainVoyager's VTC format
• SPM's Analyze or NIftI (HDR/NII) formats (either as series of 3D volumes or a single 4D file)

  fmriquality  - get some quality assurance for fMRI data
 
  FORMAT:       [fq, fmridata = ] fmriquality(images [, opts])
 
  Input fields:
 
        images      list of images (or fMRI datatype, like FMR/VTC)
        opts        optional settings
         .motcor    perform motion-correction (and re-do stats, def: false)
         .nuisreg   nuisance regressors/confounds for temp. filtering
         .pbar      1x1 xfigure::ProgressBar or xprogress object
         .prange    progress bar range (default: [0 .. 1])
         .qasheet   flag, display quality assessment sheet (default: false)
         .res       resolution (for motion detection, default: from file)
         .robfilt   flag, perform filtering robustly (default: false)
         .robmotcor flag, do motion detection robustly (default: false)
         .tempffrq  temp filter frequency cut-off (as TRs, default: 80)
         .tempfset  temp filter set, either of 'DC', {'Fourrier'}
 
  Output fields:
 
        fq          complex struct, containing masks, time courses, etc.
        fmridata    4-D data slab (motion corrected, if selected)
 
  Note: the data of one run must at least fit into memory (in single
        precision), plus some temporary arrays and, if motion correction
        is selected, with further memory allocation required!

The basic algorithm can be unpacked into the following steps:

• read in the data (memory must at least allow one run to be loaded!)
• reserve additional memory for some computation outcomes
• compute mean and standard deviation (over time) images
• detect foreground (brain) and background (black/air voxels) separately
• compute average time courses for foreground and background (global and per-slice)
• estimate the smoothness of the data (one value per volume; over time)
• create a conservative estimate of the average (temporal) standard deviation of the background
• compute a “global signal-to-noise ratio” (GlobalSNR) image (using this global background noise estimate)
• compute a “local signal-to-noise ratio” (LocalSNR) image (using each voxel's individual standard deviation)
• temporally filter timecourses (univariately) and estimate the amount of variance determined by low frequencies
• if requested, perform motion detection/realignment and re-run the computations

Additionally, at the end of steps that produce time courses, an outlier detection is performed which amasses evidence for a given volume being an outlier.

The detection of the foreground involves the following steps:

• selection of voxels for which the mean value (over time) exceeds the mean value over the entire 4D data slab
• removal of stray voxels (by one step of 3D erosion following by a back-dilation and logical AND with the original selection)
• sub-selection of the “biggest chunk” of cohesive voxels (clustering)

The detection of the background involves the following steps:

• picking the median of over-time mean values for voxels where the mean value does not exceed the mean value over the entire 4D data slab (in other words the opposite of the preliminary foreground mask!)
• selecting voxels for which the mean value is smaller than this median (very conservative background estimate)
• removal of voxels for which more than 5 per cent of values are exactly 0 (depending on the scanner type and sequence parameters as well as field homogeneity corrections, some voxels are, more or less consitently, all-0, and hence not “true background” voxels
• equally, sub-selection of the biggest chunk (to remove stray voxels)

The mentioned (conservative) estimate of the background noise is computed by

• sorting the values of the temporal standard deviation in voxels marked as background
• computing the average over the second and third quartile, so as not have voxels which are, for instance, on the fringe of the brain and, due to motion, for some portion of the run contain actual data, pollute the estimate; the idea being that the average noise picked up in air voxels should not be influenced by actual matter, even if only present in part of a run

The following criteria are being considered for the detection of outliers:

• the estimated smoothness in a given volume is further away than 6 standard deviations from the mean
• the global foreground time course is further away than 5 standard deviations from the mean
• the absolute of the 1st order derivative of the global foreground time course is further away than 5 standard deviations from the mean (detecting stark signal level shifts)
• the Mahalanobis Distance over the foreground slices' time courses is further away than 5 standard deviations from the mean (detecting patterns in time courses across slices, e.g. when partial-volume effects of motion make a volume an outlier, particularly in interleaved sequences!), this is also done with the absolute of the 1st order derivatives of the foreground slices' time courses
• the temporally fitered versions of the global foreground timecourse as well as Mahalanobis Distance over the filtered foreground slices' time courses are further away than 4 standard deviations from their respective mean

The most basic (and pre-configured) way of running fmriquality is by simply passing in the filename(s) or object of the run to check:

• using Analyze files:

  fmriquality_example1.m

  qas = fmriquality(findfiles(sessionfolder, '*.img', 'depth=1'));

• using a BrainVoyager QX FMR file:

  fmriquality_example2.m

  fmr = xff('*.fmr', 'Select FMR for which you want to check the data quality...');
  qas = fmriquality(fmr);

• also checking motion parameters with 64 TRs as the filtering cutoff, showing the result immediately:

  fmriquality_example3.m

  % selecting files with findfiles
  qafiles = fmriquality(findfiles(pwd, '*.img'));
   
  % options
  qaopts = struct( ...
      'motcor',   true, ...
      'qasheet',  true, ...
      'robfilt',  true, ...
      'tempffrq', 64);
   
  % running fmriquality
  fmriquality(qafiles, qaopts);

• create a quality assessment struct for all VTCs in an MDM, applying robust temporal filtering:

  fmriquality_example4.m

  % loop over VTCs in MDM
  for study = 1:size(mdm.XTC_RTC, 1)
   
      % perform quality assessment
      q = fmriquality(mdm.XTC_RTC{study, 1}, struct('robfilt', true));
   
      % store as VTCNAME_qasheet.mat
      save([mdm.XTC_RTC{study, 1}(1:end-4) '_qasheet.mat'], 'q');
  end

In case the QA sheet is not shown by the function (flag .qasheet set to false, which is the default), the returned struct can be visualized via a call to fmriqasheet:

fmriqasheet_example.m

% load a QA sheet (contains variable/struct q !)
load CPU4212_run3_qasheet.mat
 
% open QA sheet
fmriqasheet(q);

If the .qasheet option is not set to true, the returned variable can later be passed to fmriqasheet manually (for scripted QA-ing). This variable is of type struct and contains (at least) the following fields:

   .Dims         1x4 array, size
   .Filename     the first filename given
   .Masks        automatically detected masks (foreground, background, etc)
   .Raw          mean, stdev, and null-voxel image
   .TempFiltered re-created summary values/maps after applying temporal filtering
   .Quality      summary images trying to capture overall quality measures (SNR, CNR, etc.)
   .TC           diverse time courses

When passed to fmriqasheet, this function then creates a new figure and displays part of the information contained in the structure, which can be used to decide on whether or not a subject would likely introduce too much noise/bias at the group level.

Please be aware that the first argument must be a list of filenames; so, even for a single NIftI (NII) file, a cell array must be passed in:

fmriquality_example5.m

% using a single NII file
qas = fmriquality({'vols.nii'});

Also, please note that motion parameter estimation requires between one and several minutes (depending of the length of the run and whether or not robust motion estimation is employed). The runtime also easily increases tenfold when robust temporal filtering is used. In other words, these options are more useful and meant for situations where the fmriquality function is part of a larger script that automatically checks the quality of all incoming data instead of being used from the command line in a more interactive way!

Finally, at the moment, the “problem detection” thresholds (to mark a volume as outliers) are fixed, but are likely to change in a future version.