This file is the primary documentation for users of the fMRI Matlab software on /afs/cs/project/theo-72/fmri_core_new. The software on this directory is a self-contained core set of functions intended to be easily used on a variety of fMRI data sets. For example, the following commands in an interactive Matlab session will load the data for one subject in an fMRI experiment, display the fMRI data as an animated movie, then train a Naive Bayes classifier and test its accuracy over the training data. load data-starplus-04847-v7.mat; animateTrial(info,data,meta,3,6); %see movie [i,d,m]=transformIDM_selectTrials(info,data,meta,find([info.cond]~=0)); % seletct non-noisey trials [examples,labels,expInfo] = idmToExamples_condLabel(i,d,m); %create training data [classifier] = trainClassifier(examples,labels,'nbayes'); %train classifier [predictions] = applyClassifier(examples,classifier,'nbayes'); %test it [result,predictedLabels,trace] = summarizePredictions(predictions,classifier,'averageRank',labels); 1-result{1} % rank accuracy This documentation has the following sections: 1. An introductory tutorial 2. Description of the key data structures used by this software. 3. Documentation on user functions This software is the result of several people's efforts. The first version is due to Francisco Pereira. Other contributors include Tom Mitchell, Stefan Niculescu, Xuerui Wang, Indra Rustandi, Rebecca Hutchinson, Jay Pujara. WARNING: Some of the utilities assume your Matlab includes the statistics toolbox and the optimization toolbox. ======================================= 1. Introductory tutorial. Below is a sample session illustrating some of the most basic functions available. It is STRONGLY recommended that you type the commands in this tutorial into your Matlab window, to understand the functions and to examine the data structures and see the visualizations of the data. % load the data (assume you're using version 7 of Matlab. If you're using % version 6, load the -v6.mat file instead) load('data-starplus-04847-v7.mat') % examine the variables info, data, and meta, and read their description in % section 2 below. meta info data % draw a picture of activation for the z=4 slice of the brain for trial 4, at time 8 plotSnapshot(info,data,meta,4,4,8,0,0); % do the same as above, but for all slices of the brain plotSnapshot3D(info,data,meta,4,8,0,0); % play a movie showing the activity over time for trial 4 M=animate16Trial(info,data,meta,4); movie(M,2,.5) % play it again, twice, at .5 frames/second % plot the time course for a single voxel at coord <36,62,8> plotVoxel(info,data,meta,36,62,8) % plot the time course for just trials 3 and 4 of voxel <36,62,8> [i,d,m]=transformIDM_selectTrials(info,data,meta,[3,4]); plotVoxel(i,d,m,36,62,8); % create an info,data,meta containing just the mean activitation over the trials that require the subject to make a decision, then plot that mean activation. alltrials= 1:meta.ntrials ts= [info.actionAnswer]==0 trialsOfInterest=alltrials(ts) [i,d,m]=transformIDM_avgTrials(info,data,meta,trialsOfInterest); plotVoxel(i,d,m,36,62,8); animate16Trial(i,d,m,1); % create training examples consisting of 8 second (16 image) windows of data, labeled by whether the subject is viewing a picture or a sentence. Label 1 means they are reading a sentence, label=2 means they are viewing a picture. "examples" is a NxM array, where N is the number of training examples, and M is the number of features per training example (in this case, one feature per voxel at each of 16 time points). % collect the non-noise and non-fixation trials trials=find([info.cond]>1); [info1,data1,meta1]=transformIDM_selectTrials(info,data,meta,trials); % seperate P1st and S1st trials [infoP1,dataP1,metaP1]=transformIDM_selectTrials(info1,data1,meta1,find([info1.firstStimulus]=='P')); [infoS1,dataS1,metaS1]=transformIDM_selectTrials(info1,data1,meta1,find([info1.firstStimulus]=='S')); % seperate reading P vs S [infoP2,dataP2,metaP2]=transformIDM_selectTimewindow(infoP1,dataP1,metaP1,[1:16]); [infoP3,dataP3,metaP3]=transformIDM_selectTimewindow(infoS1,dataS1,metaS1,[17:32]); [infoS2,dataS2,metaS2]=transformIDM_selectTimewindow(infoP1,dataP1,metaP1,[17:32]); [infoS3,dataS3,metaS3]=transformIDM_selectTimewindow(infoS1,dataS1,metaS1,[1:16]); % convert to examples [examplesP2,labelsP2,exInfoP2]=idmToExamples_condLabel(infoP2,dataP2,metaP2); [examplesP3,labelsP3,exInfoP3]=idmToExamples_condLabel(infoP3,dataP3,metaP3); [examplesS2,labelsS2,exInfoS2]=idmToExamples_condLabel(infoS2,dataS2,metaS2); [examplesS3,labelsS3,exInfoS3]=idmToExamples_condLabel(infoS3,dataS3,metaS3); % combine examples and create labels. Label 'picture' 1, label 'sentence' 2. examplesP=[examplesP2;examplesP3]; examplesS=[examplesS2;examplesS3]; labelsP=ones(size(examplesP,1),1); labelsS=ones(size(examplesS,1),1)+1; examples=[examplesP;examplesS]; labels=[labelsP;labelsS]; % train a Naive Bayes classifier [classifier] = trainClassifier(examples,labels,'nbayes'); %train classifier % apply the Naive Bayes classifier to the training data (it's best to use cross validation, of course, to obtain an estimate of its true error). The returned array 'predictions' is an array where predictions(k,j) = log P(example_k | class_j). [predictions] = applyClassifier(examples,classifier); %test it % summarize the results of the above predictions. [result,predictedLabels,trace] = summarizePredictions(predictions,classifier,'averageRank',labels); 1-result{1} % rank accuracy ======================================= 2. DATA STRUCTURES see the file README-data-documentation.tex ================================================ 3. LIST OF USER FUNCTIONS This is BRIEF documentation. For more detail on function , type 'help ' in Matlab, or see the comments at the beginning of the .m file. ===================== Functions on IDM's: (IDM is short for Info,Data,Meta - see above documentation on data structures) transformIDM_selectTrials(info,data,meta,trials) Returns a copy of info,data,meta containing only the specified trials. The input parameter 'trials' is an array listing the indices of trials to be included. Example: select just trials 3 and 5 [info2,data2,meta2] = transformIDM_selectTrials(info,data,meta,[3,5]) transformIDM_selectTimewindow(info,data,meta,snapshots) Returns a copy of info,data,meta containing only the specified snapshots in time within each trial. The input parameter 'snapshots' is an array listing the indices of snapshots to be included, assuming the index of the first snapshot of each trial is 1. Example: select just time snapshots 3, 4, and 7 from each trial [info2,data2,meta2] = transformIDM_selectTimewindow(info,data,meta,[3,4,7]) transformIDM_selectVoxelSubset(info,data,meta,) Returns an IDM containing only the voxels selected in the list that is passed. Example: select voxels 3 and 5 [info2,data2,meta2] = transformIDM_selectVoxelSubset(info,data,meta,[3,5]) transformIDM_selectROIVoxels(info,data,meta,ROIlist) Returns a copy of info,data, and meta, selecting only voxels that belong to ROIs found on ROIlist. Example: [info2,data2,meta2] = transformIDM_selectROIVoxels(info,data,meta,{'CALC' 'LIT'}) returns an IDM where the data contains just the voxels belonging to CALC and LIT transformIDM_mergeTrials(info,data,meta) Transforms info, data, and meta by concatenating all trials into a single trial. The returned info has info(1).cond=-1 to indicate there is no single condition for the merged trial. It also contains info(1).snapshotCond(i) set to the condition of the trial from which the i-th snapshot was taken. The returned Meta contains meta.startsOfFormerTrials, an array listing the snapshot indices in the returned Data that correspond to the beginning snapshots of the original trials. Example: [info2,data2,meta2] = transformIDM_mergeTrials(info,data,meta) transformIDM_pairwiseAvg(info,data,meta) Returns a copy of info,data, and meta, replacing each snapshot(t) by the average of shapshot(t) and snapshot(t+1). Note this shortens the length of each trial by one snapshot. Info and data are updated accordingly. Example: [info2,data2,meta2] = transformIDM_pairwiseAvg(info,data,meta); transformIDM_selectActiveVoxact(info,data,meta,nToKeep,[conditions to consider]) Returns an IDM containing only the most active voxels across [conditions], using t-test of these conditions against condition 1 (by convention, this is 'fixation'). Conceptually, this works by taking the per condition rankings Cond 2 ... last condition being considered v(1) v(1) v(2) v(2) and -------------- sliding a window down until the number left above is the # to keep. Note that some voxels may be near the top of several of the rankings. They will get picked only once. Example: This will t-test condition 2 versus condition1, and condition 3 vs. 1, then combine selected voxels from each test [info2,data2,meta2] = transformIDM_selectActiveVoxact(info,data,meta,20,[2,3]); transformIDM_selectActiveVoxels(info,data,meta,nToKeep) Returns an IDM containing only the most active voxels from each condition. Thus, it keeps at most (#ROIS * #conditions * #nToKeep) voxels. Example: [info2,data2,meta2] = transformIDM_selectActiveVoxels(info,data,meta,20); transformIDM_avgTrials(info,data,meta,trials) Returns a copy of info,data, and meta, containing a single trial which is the average of all trials specified by the input. The input 'trials' is an array listing the indices of trials to be included in the average. If trials are of different length, prints a warning message and uses the length of the shortest trial. Example: [info2,data2,meta2] = transformIDM_avgTrials(info,data,meta,[3,6]); returns an IDM containing one trial which is the average of trials 3 and 6 transformIDM_avgTrialCondition( info,data,meta ) Returns an IDM containing one trial per condition, where the trial is the average of all the trials with that condition (with as many time points as the shortest of them). Example: [avginfo,avgdata,avgmeta] = transformIDM_avgTrialCondition(info,data,meta); transformIDM_avgVoxels(info,data,meta,superVoxels,) Returns a copy of info,data, and meta, with data defined in terms of the specified 'superVoxels'. The input superVoxels is a cell array of arrays, each of which specifies a set of voxels to be combined into one superVoxel. The ith supervoxel is assigned a pseudo-coordinate of x=i,y=1,z=1. The optional argument 'absoluteVal' determines whether to average the voxel values, or their absolute values. Set it to 1 if you wish to average the absolute values -- it defaults to 0, which specifies averaging the actual, signed values. Example: [info,data,meta] = transformIDM_avgVoxels(info,data,meta,{[3,4,5],[10,12]}); returns an IDM where the data contains just two voxels, one containing the average of voxels 3,4 and 5, and the second the avg of voxels 10 and 12. transformIDM_avgROIVoxels(info,data,meta,) Returns a copy of info,data, and meta, with each ROI replaced by a single 'superVoxel', whose activation is the mean activation of voxels in that ROI. ROIs is a cell array specifying the text names of the ROIs to include, and the sequence of the ROIs given determines the sequence of the supervoxel columns in the returned Data. If ROIs is not given, then all ROIs in meta.rois are used. The ith ROI supervoxel is assigned a pseudo-coordinate of x=i,y=1,z=1. (implemented using transformIDM_avgVoxels). **NOTE: this function assumes you have run createColToROI beforehand *** Example: [info2,data2,meta2] = transformIDM_avgROIVoxels(info,data,meta,{'LT' 'RB'}); returns an IDM where data contains just two voxels, the first giving the spatial average of voxels in LT, the second giving the mean activity of voxels in RB. transformIDM_avgVoxelSubset(info,data,meta,[voxel list,[weight list]]) Averages a selected subset of the voxels into a new IDM. The arguments are are an array with the voxel numbers to average (default all) and an array with the same size with the averaging weight (default equal for all). Example: [info2,data2,meta2] = transformIDM_avgVoxelSubset(info,data,meta,[3,4,5]); returns an IDM where the data contains one voxel, the average of voxels 3,4 and 5. transformIDM_unfold(info,data,meta) Returns an IDM where each trial gets turned into |trial| one-time point trials with the same condition. E.g. a trial with condition 3 and length 20 becomes a succession of 20 trials with condition 3 and length 1. Info,data and meta are updated accordingly. This is used in block studies to turn each time point in a block into a trial, and then each trial into an example. Example: [info2,data2,meta2] = transformIDM_unfold(info,data,meta); transformIDM_separateBlocks(info,data,meta) Returns an IDM where each trial with a given condition number acquires a different condition number. For instance, the sequence of conditions in a study with 6 conditions 2,3,2,4,3,5,6,3 would become 2,3,7,4,8,5,6,9 Info is updated accordingly. One application of this was to transform 6 into 12 categories in the data from the Categories study. Example: [info2,data2,meta2] = transformIDM_separateBlocks(info,data,meta); transformIDM_normalizeTrials(info,data,meta) Returns a copy of info,data,meta in which each trial is normalized. Normalization is done separtely for each trial, and for each voxel within each trial. Normalization consists for each voxel of calculating its mean activity over the trial, then subtracting its mean from its activity at each time point within the trial. The net effect is to reset the mean to zero, while preserving the signal variation of the voxel over time during the trial. May be useful for comparing trials, because it can remove the effect of long-term drift in signal magnitude. Example: [info2,data2,meta2] = transformIDM_normalizeTrials(info,data,meta) transformIDM_smoothingKR(info,data,meta) Smooths the IDM temporally by performing Kernel Regression with gaussian kernels, on a trial by trial basis. WARNING: for many voxels, the residuals after smoothing (data - smoothed data) are strongly correlated (i.e. they do not look like a sample from a mean 0 gaussian). This leads to a few voxels being undersmoothed, as fitting the noise in points around the one being left out will lead to a good fit there as well. It's not a major issue, in that most voxels seem appropriately smooth, and the worse that can happen is that a voxel will be undersmoothed, rather than oversmoothed (which could destroy information). Example: [info2,data2,meta2] = transformIDM_smoothingKR(info,data,meta,[conditions]) where the optional argument is a list of conditions to smooth (defaults to all). ===================== Functions on Examples: idmToExamples_fixation(info,data,meta) Produces array of examples from the info,data,meta, labeling examples as either the fixation condition or some other condition. Example: [examples,labels,expInfo] = idmToExamples_fixation(info,data,meta); idmToExamples_condLabel(info,data,meta) Convert info, data, meta structure to example, lables structure, labeling examples as the condition numbers from info. Exampels: [examples,labels,expInfo] = idmToExamples_condLabel(info,data,meta); mergeExamples(examples1, labels1, expInfo1, examples2, labels2, expInfo2) Merge two example sets into a large one, and update the experiment information. Example: [examples, labels, expInfo] = mergeExamples(examples1, labels1, expInfo1, examples2, labels2, expInfo2) ===================== Functions on Classifiers: trainClassifier(Examples,Labels,classifierType) Trains on a set of examples, producing as output a classifier. For now, use this only to train a NaiveBayes classifier, by using the value 'nbayes' for classifierType. Example: classifier = trainClassifier(examples, labels, 'nbayes'); applyClassifier(Examples,classifier) Use a learned classifier (from trainClassifier) to label new examples. For nbayes classifier, outputs an array with one row per example, and one column per class. The columns are ordered in increasing order with the class label. The value of predictions(k,j) is log P(class_j | example_k). Example: [predictions] = applyClassifier(examples,classifier); summarizePredictions(predictions,classifier,'averageRank', truelabels) Given a set of examples, a classifier returned by trainClassifier, the predictions returned by applyClassifier, and a column vector of truelabels, returns a summary of the predictions. The most useful part of this is predictedLabels. The value of predictedLabels(k) is the label predicted for example k. Example: [result,predictedLabels,trace] = summarizePredictions(scores,trainedClassifier,'averageRank',testLabels) ===================== Displaying data: HINT: new plots reuse the existing Matlab plot window. To preserve existing plots and create a new one, give the Matlab command 'figure' before calling the function. plotVoxel(i,d,m,x,y,z) Display the entire time series for voxel , and the series of trial conditions, concatenating all trials in i,d,m. The voxel activation is plotted in blue. In addition, the condition number of the trial (info.cond) is plotted in red. Example: plotVoxel(i,d,m,46,54,10) % plot timecourse for voxel <46,54,10> and conditions plotSnapshot(i,d,m,trialNum,z,n,minActiv,maxActiv) Plot an image showing slice z of the n-th snapshot of trial numbered trialNum. If you want to pass in the maximum and minimum activities (i.e., the activity levels that will display as color intensities 1 and 128), then set maxActiv, minActiv accordingly. If you set these both to zero, they will be computed from the data. Example: display slice z=3 of snapshot 4 of trial 2 plotSnapshot(i,d,m,2,3,4,0,0) plotSnapshot3D(i,d,m,trialNum,n,minActiv,maxActiv) Same as plotSnapshot, but displays up to 16 z slice images at once. Z slices are laid out by z values in the following arrangement: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Example: display first snapshot of second trial. img = plotSnapshot3D(i,d,m,2,1,0,0); plotTrial(i,d,m,trialNum,) Plot the voxel time courses for trial number trialNum, one z slice at a time, for the given i,d,m. Lays out voxel plots showing their geometric relationship. Example: viewTrial(i,d,m,12) % plots time courses for trial number 12 animateTrial(i,d,m,trialNum,z) Create a movie of brain activity for one z slice of trial trialNum of the IDM, and display it. Returns the movie matrix. Example: Show and return a movie for trial 3, brain slice z=7 M=animateTrial(i,d,m,3,7); movie(M,3,2); % replay it 3 times, at 2 frames/second animate6Trial(i,d,m,trialNum,z) Just like animateTrial, but simultaneously display six z slices (z through z+5) Example: Show and return a movie for trial 3, brain slices z= 9 - 14 M=animate6Trial(i,d,m,3,9); movie(M,3,2); % replay it 3 times, at 2 frames/second animate16Trial(i,d,m,trialNum) Just like animateTrial, but simultaneously display all 16 z slices Z slices are laid out in the following format: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Example: Show and return a movie for trial 3, all brain slices M=animate16Trial(i,d,m,3); movie(M,3,2); % replay it 3 times, at 2 frames/second ===================== Utility functions: getVoxelTimecourse(info,data,meta,x,y,z) Returns a column vector containing the entire time course of voxel throughout the IDM getConditionTimecourse(info,data,meta) Returns a column vector containing the condition numbers at each time instant throughout the IDM getROI(meta,roiNameString) Returns the struct from meta.rois that corresponds to 'roiNameString' Example: getROI(meta,'LT')