From your terminal type
echo $FSLDIRto find out what the FSL path is, and then
addpath('[whatever the FSL path is]/etc/matlab');from MATLAB. To load the data and mask files
[data4D] = read_avw('dataset_small.nii.gz'); [mask3D] = read_avw('dataset_small_atlas.nii.gz');Note of caution:
Please note that these functions do not read or write orientation information in the NIfTI files (i.e. qform and sform) and so can swap left and right for the specific NIfTI files that you are using if you are not careful; you should be alert to this, and examine the starting files (with fslhd and fslview) and/or plot a few slices from the mask3D where there should be asymmetries to see if it has occurred. You can deal with this within MATLAB and there are FSL functions you can use to modify the NIfTI files a priori (fslcpgeom/fslorient); for more on the latter, please read Fslutils and Orientation Explained.
You can get it from
http://sourceforge.net/projects/niftilib/files/or simply use the version bundled with this toolbox. After unpacking you should
addpath('[whatever path to unpacked folder]/matlab');from MATLAB. To load the data and mask files
system('gunzip dataset_small.nii.gz'); tmp = nifti('dataset_small.nii'); data4D = tmp.dat(:,:,:,:); system('gzip dataset_small.nii'); system('gunzip dataset_small_atlas.nii.gz'); tmp = nifti('dataset_small_atlas.nii'); mask3D = tmp.dat(:,:,:); system('gzip dataset_small_atlas.nii');The calls to gzip and gunzip are necessary because the package does not handle gzippped NIfTI files (if your data files end in .nii these calls will not be necessary).
size(data4D) size(mask3D)to see their dimensions. The assumption for this tutorial is that the 160 volumes in the data4D matrix correspond to 160 trials, e.g. each volume might be the average image around the peak of the hemodynamic response in a slow event related design, or a volume of beta coefficients deconvolved from a faster design.
The next step is to make labels for each of those volumes. We can load regressor files directly
regressorsCondition = load('dataset_small_regressors_condition.txt'); regressorsRun = load('dataset_small_regressors_run.txt');(which are both 160x4, as there are 4 conditions and 4 runs), and convert them into 160 entry vectors
labels = sum(repmat([1 2 3 4],160,1) .* regressorsCondition,2); labelsRun = sum(repmat([1 2 3 4],160,1) .* regressorsRun,2);
meta = createMetaFromMask(mask3D,'radius',1);using a radius of 1 (corresponding to a 3x3x3 searchlight), and the result is a structure containing that meta-information (described in the next section if you haven't used it before).
We then use the meta structure to generate a set of examples, each of which is a vector containing all the information in the brain mask voxels only for a particular volume. This saves a lot of space and also helps speed up calculations.
examples = createExamplesFrom4D(data4D,meta); size(examples)Note that examples is a 160x896 array, as there are only 896 voxels in the brain mask used.
By now you've seen data formats such as nifti that store entire 3D volumes over time. This is somewhat wasteful, as we are often interested in a fraction of those voxels (e.g. those in a mask covering cortex). The CMU format stores data as an array with #examples x #voxels, including only those mask voxels.
The best way of introducing this is with an example. Let's say we a 4D dataset with a 3D volume for each of 10 time points, and the 3D volume is 4x4x4
dataset = randn(4,4,4,10);Even though these are little volumes, we would like to look at the tiny 2x2x2 brain inside, as defined by this mask
mask = zeros(4,4,4); mask(2:3,2:3,2:3) = 1;If you want to see what it looks like, you can disp(mask) or
clf; for i = 1:4 subplot(1,4,i); imagesc(mask(:,:,i),[0 1]); axis square; endThe "meta" structure is produced for a mask over the 3D volume that we would like to use, using
tinymeta = createMetaFromMask(mask);If you type
tinymetayou'll see various fields. For now, please look at
dimensions - the dimensions of the imaging volume dimx,dimy,dimz - same, if you quickly want to get each one indicesIn3D - m element vector these are just the linear indices into a [dimx] x [dimy] x [dimz] matrix of the value 1 voxels in the mask.If you type
mask(tinymeta.indicesIn3D)you should thus get a vector with all 1s. This can also be used to pick the values of mask 1-voxels in any volume in the dataset, e.g.
volume = dataset(:,:,:,1); values = volume(tinymeta.indicesIn3D);
tinyexamples(1,:)and make sure the appropriate voxels in dataset(:,:,:,1) are picked up.
EXERCISE 1 SOLUTION: nvoxels = length(tinymeta.indicesIn3D); tinyexamples = zeros(10,nvoxels); for i = 1:10 volume = dataset(:,:,:,i); tinyexamples(i,:) = volume(tinymeta.indicesIn3D); end
colToCoord - #voxels x 3 matrix vector = coordToCol(i,:) places the 3D coordinates of voxel i in vector coordToCol -
x x matrix i = colToCoord(x,y,z) gives the column of voxel with 3D coordinates x, y and z If i = 0 then the voxel is not in the data matrix.
EXERCISE 2 SOLUTION voxel = tinymeta.coordToCol(2,2,2); values = tinyexamples(:,voxel) coordinates = tinymeta.colToCoord(8,:)
numberOfNeighbours - #voxels x 1 numberOfNeighbours(i) is the number of neighbours a voxel has voxelsToNeighbours - #voxels x (2*radius+1)^3 neighbours = voxelsToNeighbours(i,1:numberOfNeighbours(i)); are the neighbours of voxel i (not including i). Note that row i in voxelsToNeighbours will be filled with numbers only if voxel i has neighbours in every possible direction. If not, the entries above numberOfNeighbours(i) are just NaN.
EXERCISE 3 SOLUTION voxel = tinymeta.coordToCol(2,2,2); number = tinymeta.numberOfNeighbours(voxel); neighbours = tinymeta.voxelsToNeighbours(voxel,1:number)
tinymeta2 = createMetaFromMask(mask,'radius',2);This is very similar to tinymeta, except for one of the fields. Which one?
EXERCISE 4 SOLUTION size(tinymeta2.voxelsToNeighbours) will show you that, in a searchlight with radius 2 (5x5x5 cube) you have up to 124 neighbours.
load('dataset_small.mat'); metaYou will notice that, in addition to the fields discussed above, this meta has 3 more fields
nrois roiIDs roiColumns
This was obtained by giving createMetaFromMask a mask volume containing not just 1s where the brain was present and 0s everywhere else but, instead, a numeric ROI ID in each voxel corresponding to the ROI that voxel belonged to. In the tutorial dataset there are 4 ROIs, with IDs stored in meta.roiIDs, and meta.roiColumns contains, for each of those 4, the column numbers in examples containing the voxels belonging to that ROI. We can use these numbers to find their respective 3D coordinates, by using them as indices into meta.colToCoord. Absent ROI IDs in a 3D mask, meta will contain a single ROI.