Voxel and image spaces in CANlab fmri_data objects

This brief guide looks at how images are stored as matrices of different dimensions in fmri_data objects. It covers:
Table of Contents
When we load a set of images from .nii files, it automatically applies a default "whole brain mask" to save space.
The rationale is that there's no point in saving empty values outside the brain.
% Load a pre-canned example dataset
imgs = load_image_set('emotionreg', 'noverbose');
Defaults settings have been modified by file(s): /Users/f003vz1/Dropbox (Dartmouth College)/Matlab_code_external/spm12/spm_my_defaults.m Modified fields: stats
sampleto = '/Users/f003vz1/Documents/GitHub/CanlabCore/CanlabCore/Sample_datasets/Wager_et_al_2008_Neuron_EmotionReg/Wager_2008_emo_reg_vs_look_neg_contrast_images.nii'
Defaults settings have been modified by file(s): /Users/f003vz1/Dropbox (Dartmouth College)/Matlab_code_external/spm12/spm_my_defaults.m Modified fields: stats
This is what the default mask looks like (all the colored voxels are in-mask).
Only values within this mask will be stored in the object.
figure; axis off; montage(imgs.mask, 'trans');
Setting up fmridisplay objects sagittal montage: 4895 voxels displayed, 269358 not displayed on these slices sagittal montage: 4728 voxels displayed, 269525 not displayed on these slices sagittal montage: 4794 voxels displayed, 269459 not displayed on these slices axial montage: 38512 voxels displayed, 235741 not displayed on these slices axial montage: 41445 voxels displayed, 232808 not displayed on these slices
The actual image data are stored in the object in the .dat field, which is a matrix of [voxels x images] dimension. Here, we have [31686 in-mask voxels x 30 images] stored in the object:
size(imgs.dat)
ans = 1×2
352328 30

Voxel sets associated with the image object

There are three sets of voxels associated with the image object:
  1. The entire image set of E voxels
  2. The "in-mask" set of M voxels (excluding voxels outside the mask used when creating the image)
  3. The "data list" set of V voxels (excluding voxels with missing/empty values (0 or NaN) for all images
The rationale is that we don't need to save values for thousands of voxels if they are empty in all images. (0 or NaN, zero is coded as missing/empty).
The method remove_empty( ) finds and removes voxels empty in all images, and stores a logical indicator for which were removed in the .removed_voxels field. Before running this, the .dat field has a data matrix of [M x number_of_images]. After running, the .dat field now has a data matrix of [V x number_of_images].The .removed_voxels field is [M x 1].
The method replace_empty( ) uses the the .removed_voxels field to reconstruct the .dat field to the original [M x number_of_images]. Empty voxels are filled in with zeros in all images.
Sometimes, it's advantageous to work only with the "data list", where the data matrix has at least one valid value for each voxel. Other times, as when you're matching up images to apply a mask or concatenate data, you need the full set of M voxels in the original mask.
Finally, other times you may need a dataset with one value for the entire set of E voxels in the image. It's not possible to recover values outside the mask, as they were omitted during object creation. However, you can use the method reconstruct_image to reconstruct the entire 3-D (for an object with one image) or 4-D (for an object with multiple images) data matrix. This fills empty values with zero. Then, you can vectorize this 3-D or 4-D matrix.

The .volInfo field stores voxel and image space information

The .volInfo field has information aboutt the mask, space, voxel sizes, origin, and transformations from "image" (voxel) to "world" (mm) space.
imgs.volInfo
ans = struct with fields:
fname: 'REMOVED: CHANGED SPACE' dim: [47 56 31] dt: [16 1] pinfo: [3×1 double] mat: [4×4 double] n: [1 1] descrip: 'Space of /Users/f003vz1/Documents/GitHub/CanlabCore/CanlabCore/Sample_datasets/Wager_et_al_2008_Neuron_EmotionReg/Wager_2008_emo_reg_vs_look_neg_contrast_images.nii' private: [1×1 nifti] nvox: 81592 image_indx: [81592×1 logical] wh_inmask: [35676×1 double] n_inmask: 35676 xyzlist: [35676×3 double] cluster: [35676×1 double]
imgs.volInfo.mat contains the matrix with the origin (zero point in mm), voxel sizes, and any rigid-body rotations applied to the image. This is the same format as used by SPM and stored in SPM ".mat" files.
imgs.volInfo.mat
ans = 4×4
-3.4375 0 0 82.5000 0 3.4375 0 -116.8750 0 0 4.5000 -54.0000 0 0 0 1.0000
You shouldn't have to worry about these much, as the object methods should handle them "under the hood". But sometimes it's valuable to know.
imgs.volInfo.dim contains the 3-D dimensions of the image in voxels, in x (left-to-right), y (back-to-front), and z (bottom-to-top).
imgs.volInfo.dim
ans = 1×3
47 56 31
imgs.volInfo.xyzlist contains a list of the voxel coordinates [x y z] for each in-mask voxel (an [M x 3] matrix)
and imgs.volInfo.cluster contains an [M x 1] integer vector, where the integers represent contiguous "blobs" estimated with spm_cluster.m.
These are used in various operations, including the transformation of an fmri_data object into a region-class object (with region( )),which stores data grouped by contiguous region.
imgs.volInfo.xyzlist(1:5, :)
ans = 5×3
18 9 1 19 9 1 28 9 1 29 9 1 30 9 1
imgs.volInfo.cluster(1:5, :)
ans = 5×1
1 1 1 1 1

An example looking at the number of voxels in each voxel list

Let's take a look at our example:
  1. Here there are 81,592 voxels in the entire image (E), including "out of mask" voxels
  2. There are 35,676 "in-mask" voxels (M)
  3. We start off with 35,676 voxels in the "data list"
% 81592 in the entire image, including "out of mask" voxels
imgs.volInfo.nvox
ans = 81592
% 35676 in the brain mask (with possible data, but some voxels may be empty)
imgs.volInfo.n_inmask
ans = 35676
% 35676 voxels in the data list in .dat
size(imgs.dat, 1)
ans = 35676
Now we'll use remove_empty to exclude empty/missing voxels from the data list.
We still have 35,676 voxels in the "data list" because the images contain values for all these voxels.
% 35676 voxels with non-empty valid data (non-zero, non-nan) for at least 1 image
imgs = remove_empty(imgs);
size(imgs.dat, 1)
ans = 35676

Masking

Now we'll apply a gray-matter mask, excluding voxels in the ventricles and gray matter, which will be removed from the image object. These values are discarded, so we can't recover them later.
imgs = apply_mask(imgs, which('gray_matter_mask.nii'));
Defaults settings have been modified by file(s): /Users/f003vz1/Dropbox (Dartmouth College)/Matlab_code_external/spm12/spm_my_defaults.m Modified fields: stats
figure; axis off; montage(mean(imgs), 'trans');
Setting up fmridisplay objects sagittal montage: 1103 voxels displayed, 30583 not displayed on these slices sagittal montage: 1061 voxels displayed, 30625 not displayed on these slices sagittal montage: 1010 voxels displayed, 30676 not displayed on these slices axial montage: 10380 voxels displayed, 21306 not displayed on these slices axial montage: 10986 voxels displayed, 20700 not displayed on these slices
Now, we have a reduced number of voxels with valid data, as we've masked some out:
size(imgs.dat) % now we have only 31,686 voxels with data because we have masked some out
ans = 1×2
31686 30
length(imgs.removed_voxels) % this is still the length of the original mask
ans = 35676
The sum of the removed_voxels logical vector reflects the number removed. The number removed and the number remaining must always sum to M (number in mask) for any valid object:
sum(imgs.removed_voxels) % this is how many voxels are excluded from current imgs.dat
ans = 3990
size(imgs.dat, 1) + sum(imgs.removed_voxels) % current voxels + removed must = total number of voxels in mask
ans = 35676

Resampling

Resampling uses the volInfo.mat field to calculate the rotations and interpolation (i.e., changes in voxel size) needed to sample an image object into the same space as a reference image object, so that all the voxels line up (i.e., are located at the same brain coordinates) across objects.
This is essential when applying masks, patterns, and atlas regions of interest, and when concatenating data for analysis of combined datasets.
Resampling is often done "under the hood" when you use various methods, like apply_mask( ) or cat( ).
Let's load an atlas-class object with a number of regions of interest (489 ROIs) and resample our dataset to that atlas, so the voxels line up with the atlas labels:
atl = load_atlas('canlab2018_2mm');
Loading atlas: CANlab_combined_atlas_object_2018_2mm.mat
imgs = resample_space(imgs, atl);
imgs.volInfo.nvox
ans = 902629
size(imgs.dat)
ans = 1×2
352328 30
Now there are 902,629 voxels in the entire image (E) and 352,328 voxels in the mask (M). This matches the list in the atlas object.

Extracting data from matched objects

Since the voxels match in location across the dataset object (imgs) and atlas (atl), we can identify voxels in any given region of the atlas, and simply grab those from the .dat field in the imgs object.
% select an arbitrary region from the atlas, #23, which is area 55b
wh_vox = atl.dat == 23;
atl.labels{23}
ans = 'Ctx_55b_L'
% get the [voxels x images] matrix of data from our test dataset, for
% voxels in area 55b
roi_data = imgs.dat(wh_vox, :);
whos roi_data
Name Size Bytes Class Attributes roi_data 340x30 81600 double
There are other ways of extracting data as well. This example is just to illustrate voxel lists and spaces.
Some other methods include:
select_atlas_subset(atl) to get one or more regions by name or number from the atlas object
extract_roi_averages(imgs, atl) to get ROI averages or pattern expression data from all the regions in the atlas.
data_table = extract_data(atl, imgs) to do the same (this is an atlas-object method)
extract_roi_averages(atl, imgs) to do the same (this is an atlas-object method)
DAT = extract_measures_batch(imgs); to extract a bunch of signature pattern responses and ROI means from several parcellations.