Writing Analyses

mrLoadRet is designed to make it relatively easy to write an entirely new data analysis stream. The quickest way to do this is to script it and use mrDispOverlay to display the results. A complete implementation that is fully integrated with the mrLoadRet GUI is only a bit more complicated. For a couple of relatively simple examples, see:

  • mrTools/mrLoadRet/Analysis/timeSeriesStats.m
  • mrTools/mrLoadRet/Analysis/AverageTSeries

There are typically two kinds of data analyses. The first takes one or more time series as input and produces one or more time series as output. Examples of this include motion compensation, concatenation, and averaging time series. The second type takes one or more time series as input and produces one or more maps (statistical parameter maps or overlays) as output. Examples of this include time series statistics, corAnal for analyzing retinotopy data, and GLM for analyzing event-related data. Each analysis is typically controlled by a handful of parameters.

For each new analysis, you are free to parameterize it in any way you wish. However, it will be easier for you if you adopt a particular convention for the parameters structure as outlined below. Typically, we provide a GUI to allow users to choose values for the parameters. You are free to write that GUI however you wish. For some of our analyses, we used Matlab's GUIDE to create the GUI. For other analyses, we used mrParamsDialog which automatically generates a GUI.

There are, however, some constraints on how you write new analysis code. First, you are strongly encouraged to follow our programming conventions, for example, using viewGet and viewSet to access and modify mrLoadRet's data structures. Second, you will have to follow the conventions for adding time series outputs and/or maps/overlays to the mrLoadRet data structures. Third, you have to provide mrLoadRet with the information it needs to reconcile inconsistencies in the analysis parameters. This is discussed in detail below but the basic issue is as follows. Let's say that you perform an analysis like timeSeriesStats on all the scans in a particular group. Then later (perhaps weeks later) you delete one of the scans from that group and add new scan to that group. What should we do with the previously computed (and saved analysis)? We could dump it and make you recompute it but only 2 of perhaps many scans have been changed. We implemented mrLoadRet such that it notices that some scans have changed so you handle such inconsistencies efficiently. The code that we have written for this is a bit complicated because it interacts with how the analysis parameters are specified and how the results of the analysis are saved, but it will not be difficult for you if you follow our conventions.

Saving results

The last few lines of averageTSeries.m provide an example of how to save the results of an analysis that produces time series as the output. The variable 'aveTSeries' is a 4D array (a time series of volumes), and 'view' is the mrLoadRet view from which the average was computed. The variable 'baseScan' is a string specifying one of the scans that went into the average. The variable 'nFrames' is the number of frames in the time series, 'description' is a string, 'vol2mag' and 'vol2tal' are 4×4 transform matrices. The critical line of code is the last one, calling 'saveNewTSeries' which not only saves it as a nifti file but also inserts information about the new time series into the mrLoadRet data structures.

viewAverage = newView(viewGet(view,'viewType'));
 
% Save aveTSeries (using header of 1st scan on scanList as the template for
% the nifti header), and add it as a new scan.
scanParams.fileName = [];
scanParams.junkFrames = 0;
scanParams.nFrames = nFrames;
scanParams.description = description;
scanParams.vol2mag = vol2mag;
scanParams.vol2tal = vol2tal;
hdr = cbiReadNiftiHeader(viewGet(view,'tseriesPath',baseScan));
[viewAverage,tseriesFileName] = saveNewTSeries(viewAverage,aveTSeries,scanParams,hdr);

The last few lines of timeSeriesStats.m provide an example of how to add an analysis structure, and save it, for an analysis that produces maps as the output. The analysis computes several maps (mean of each voxel over time, median, standard deviation, max frame-to-frame difference, and max difference between each frame and the median). Each of these overlays is itself a data structure. See timeSeriesStats.m for details.

% Make analysis structure
tsStats.name = 'timeSeriesStats';  % This can be reset by editAnalysisGUI
tsStats.type = 'timeSeriesStats';
tsStats.groupName = params.groupName;  % Specified earlier in timeSeriesStats.m
tsStats.function = 'timeSeriesStats';
tsStats.guiFunction = 'timeSeriesStatsGUI';  % Name of function for the GUI
tsStats.params = params;  % Parameters structure
 
% Install the analysis in the view along with its overlays
view = viewSet(view,'newanalysis',tsStats);
view = viewSet(view,'newoverlay',tsMean);
view = viewSet(view,'newoverlay',tsMedian);
view = viewSet(view,'newoverlay',tsStd);
view = viewSet(view,'newoverlay',tsMaxFrameDiff);
view = viewSet(view,'newoverlay',tsMaxMedianDiff);
 
% Save it
saveAnalysis(view,tsStats.name);

Reconcile functions

Every analysis must have a reconcile function. The main purpose of the reconcile function is to match time series file names with scan numbers. For an analysis that produces maps/overlays as the output, you can use as examples either corAnal or timeSeriesStats. The corAnal analysis uses 'corAnalReconcileParams.m' and the timeSeriesStats analysis adopts a particular convention for its parameters structure so it uses 'defaultReconcileParams.m'. For an analysis that produces time series as the output, you can use as examples either 'averageTSeriesReconcileParams.m' or 'concatTSeriesReconcileParams.m'.

We'll start with timeSeriesStatistics because it is the simplest of these examples. The following script will run the timeSeriesAnalysis on a 'view' structure with group of scans called 'myGroup':

params.groupName = 'myGroup';
params.scanList = [1,4];
view = timeSeriesStats(view,params);

Note that other than the groupName, the only parameter for this analysis is a list of the scans upon which to run it. The resulting analysis structure looks like this:

>> view.analyses{1}
           curOverlay: 5
                    d: []
                 date: '12-Jul-2009 14:02:34'
             function: 'timeSeriesStats'
            groupName: 'myGroup'
          guiFunction: 'timeSeriesStatsGUI'
        mergeFunction: 'defaultMergeParams'
                 name: 'timeSeriesStats'
             overlays: [1x5 struct]
               params: [1x1 struct]
    reconcileFunction: 'defaultReconcileParams'
                 type: 'timeSeriesStats'

and the params sub-structure now looks like this after running it:

>> view.analyses{1}.params
      groupName: 'myGroup'
       scanList: [1 4]
    tseriesFile: {'090522181817.nii'  '090522184252.nii'}

Notice that the resulting params structure now has a 'tseriesFile' field in addition to the 'groupName' and 'scanList' fields that were set before running the analysis. This is how mrLoadRet keeps track of which time series files were included in the analysis. At the time that this analysis was done, scan number 1 in the 'myGroup' group corresponded to '090522181817.nii' and scan number 4 corresponded to '090522184252.nii'. But let's say that you delete scan 1:

view = deleteScans(view,[1],'myGroup');

When you do this, the correspondence between scan numbers and time series files is automatically updated:

>> view.analyses{1}.params
      groupName: 'myGroup'
       scanList: 3
    tseriesFile: {'090522184252.nii'}

Specifically, what happened here is that 'deleteScans' called 'viewSet(view,'deleteScan',scanNum)' which in turn called 'reconcileAnalyses(view)'. The 'reconcileAnalyses' function is called also by 'viewSet(view,'newScan',fileName,description)', and it loops through all of the analyses in all of the views and evaluates the reconcile function ('defaultReconcileParams' in this case) to fix the correspondence between scan numbers and tseries files. In the process, if there are any other parameters (which there aren't in this case) that vary from scan-to-scan, then it updates those as well. The meat of reconcile function is to generate a list of current time series filenames associated with a group, then to loop through those time series filenames and look for them in 'params.tseriesFile'. For each one that is found, the corresponding parameters are copied into a new params structure. We also use the reconcile parameters function either to generate a default parameters structure or to ensure that it has the required fields.

The reconcile function also reconciles the results of the analysis (the maps/overlays) along with the analysis parameters. In the preceding example, let's look at one of the overlay structures before and after deleting the scan:

>> view.analyses{1}.overlays(1).data
   data: {[64x64x22 double] []  []  [64x64x22 double]  []  []  []  []}
>> v = deleteScans(view,[1],'myGroup');
>> view.analyses{1}.overlays(1).data
   data: {[]  []  [64x64x22 double]  []  []  []  []}

Here's what all of that looks like for 'corAnalReconcileParams':

function [newparams,newdata] = corAnalReconcileParams(groupName,params,data)
% params = corAnalReconcileParams(groupName,[params],[data])
%
% Checks for consistency between corAnal parameters and current status of
% group (nscans and tseries filenames). Reconciles params by matching
% tseries filenames. Also reorders data arrays along with params.
%
% params: Optional initial parameters (see corAnal).
% default if not passed: initialize new params with default values.
%
% data: cell array of data arrays
% default if not passed: then returns cell array of empty matrices
%
% djh 7/2006
 
% Note that you can get groupNum and nScans with [] as the first argument
% to viewGet, i.e., without specifying a view.
groupNum = viewGet([],'groupNum',groupName);
nScans = viewGet([],'nscans',groupNum);
 
% Initialize newparams
newparams.groupName = groupName;
newparams.recompute = zeros(1,nScans);
newparams.ncycles = zeros(1,nScans);
newparams.detrend = cell(1,nScans);
newparams.spatialnorm = cell(1,nScans);
newparams.tseriesfile = cell(1,nScans);
 
% Choose default values for the parameters
for scan = 1:nScans
  newparams.detrend{scan} = 'Highpass';
  newparams.spatialnorm{scan} = 'Divide by mean';
  tseriesfile = viewGet([],'tseriesFile',scan,groupNum);
  newparams.tseriesfile{scan} = tseriesfile;
  newdata{scan} = [];
end
 
% Initialize newdata (and data as well if it doesn't already exist)
if ieNotDefined('data')
  data = cell(1,nScans);
end
newdata = cell(1,nScans);
 
% Find scans with tseries files that match those specified in
% params.tseriesfile. Use only those scans and the corresponding params.
% params.tseriesfile is a cell array of strings specifying tseries
% filenames. Or it can be 'any' to allow any match with the tseries files.
% The 'any' option is useful so that you can run the analysis on multiple data 
% sets using the same params.
if ~ieNotDefined('params')
  for scan = 1:nScans
    tseriesfile = viewGet([],'tseriesFile',scan,groupNum);
    if strcmp(params.tseriesfile,'any')
      match = scan;
    else
      match = find(strcmp(tseriesfile,{params.tseriesfile{:}}));
    end
    if match
      newparams.recompute(scan) = params.recompute(match);
      newparams.ncycles(scan) = params.ncycles(match);
      newparams.detrend{scan} = params.detrend{match};
      newparams.spatialnorm{scan} = params.spatialnorm{match};
      newparams.tseriesfile{scan} = tseriesfile;
      if (length(data) >= match)
          newdata{scan} = data{match};
      end
    end
  end
end

Merge functions

Any analysis that produces maps/overlays as the output (i.e., if it constructs an analysis structure) must also have a merge function. The main purpose of the merge function is to figure out what to do when saving an analysis (or when saving an overlay) if a previously saved version of the analysis (or overlay) already exists. As examples, the corAnal analysis uses 'corAnalMergeParams.m' and the timeSeriesStats analysis uses 'defaultMergeParams.m'.

Here's what it looks like for 'corAnalMergeParams'. Note that the 'recompute' parameter determines whether or not to perform the analysis on each scan.

function [mergedParams,mergedData] = corAnalMergeParams(groupName,oldParams,newParams,oldData,newData)
%
% [mergedParams,mergedData] = corAnalMergeParams(groupName,oldParams,newParams,[oldData],[newData])
%
% Merges corAnal params. Called by saveOverlay or saveAnalysis when file already 
% exists and 'merge' option is selected. Calls corAnalReconcileParams to make sure 
% that the merged params and data are consistent with the tseries files.
%
% oldParams: corAnal params structure in saved file.
% newParmas: corAnal params structure corresponding to newly computed analysis.
% oldData: cell array of data arrays corresponding to oldParams.
%            Default: []
% newData: cell array of data arrays corresponding to newParams.
%            Default: []
%
% mergedParams: uses newParams if recompute param is non-zero. Otherwise, uses oldParams.
% mergedData: corresponding merged data.
%
% djh 5/2007
 
if ieNotDefined('oldData')
  oldData = [];
end
if ieNotDefined('newData')
  newData = [];
end
 
% Initialize mergedParams and mergedData
mergedParams = newParams;
mergedData = newData;
 
% Get nScans
groupNum = viewGet([],'groupNum',groupName);
nScans = viewGet([],'nscans',groupNum);
 
% Loop through scans, copying in the oldParams and oldData as long as it wasn't
% just recomputed.
for scan = 1:nScans
  if ~newParams.recompute(scan)
    mergedParams.recompute(scan) = oldParams.recompute(scan);
    mergedParams.ncycles(scan) = oldParams.ncycles(scan);
    mergedParams.detrend(scan) = oldParams.detrend(scan);
    mergedParams.spatialnorm(scan) = oldParams.spatialnorm(scan);
    mergedParams.tseriesfile(scan) = oldParams.tseriesfile(scan);
    if iscell(oldData)
      mergedData{scan} = oldData{scan};
    end
  end
end

Default Reconcile and Merge Param functions

The above describes how you write your own reconcile and merge params functions, but sometimes you may just need a bare bones functions that makes sure that your analysis structure keeps track of any scans that are deleted or added so that your saved analysis will always correctly match the scans in your group. To do that, you can just specify defaultReconcileParams and defaultMergeParams as your reconcile and merge functions. An example can be found either in concatTSeries or eventRelated:

erAnal.name = params.saveName;
erAnal.type = 'erAnal';
erAnal.groupName = params.groupName;
erAnal.function = 'eventRelated';
erAnal.reconcileFunction = 'defaultReconcileParams';
erAnal.mergeFunction = 'defaultMergeParams';
erAnal.guiFunction = 'eventRelatedGUI';
erAnal.params = params;
erAnal.overlays = r2;
erAnal.curOverlay = 1;
erAnal.date = dateString;
view = viewSet(view,'newAnalysis',erAnal);

Note that your params, should have the following field:

params.scanList

which contains an array with which scan numbers you have done the analysis over. The defaultReconcileParams function will fill in field that contains the corresponding tseriesFiles which are used for matching the analysis with the current scans (tseriesFilenames are unique, but scan numbers can get changed as you delete and insert scans).

If you have created your params using the function mrParamsdialog then defaultReconcileParams will do some extra validation against the paramsInfo field of params.