%% 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:
% The basics of how to understand the masking procedure, data matrices, and voxel lists in fmri_data objects
% Where space information mapping from matrix to 3-D brain space is stored
% How to reconstruct an image into 3-D or 4-D matrix, resample one object to the space of another, and extract data from regions of interest specified by one object (e.g., an atlas/parcellation) from another (a test image dataset).
% Table of Contents
% Voxel sets associated with the image object
% The .volInfo field stores voxel and image space information
% An example looking at the number of voxels in each voxel list
% Masking
% Resampling
% Extracting data from matched objects
% 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');

%
% 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');

% 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)  

%
%% Voxel sets associated with the image object
% There are three sets of voxels associated with the image object:
% The entire image set of E voxels
% The "in-mask" set of M voxels (excluding voxels outside the mask used when creating the image)
% 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

% 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

% Here, the voxel size is 3.4375 and the origin is [82.5000 -116.8750  -54.0000].
% The negative sign of the first element (-3.4375) indicates that the images were stored in "radiological format" and flipped left-right when loading them from disk to store them in "neurological format" (right side of image is right side of brain) in the object.
% The zeros in the off-diagonals indicate that there were no rotations applied when loading the images from .nii/.img
% The bottom right element (4, 4) is the SPM scaling factor, also applied when loading the images from .nii/.img
% 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

% 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, :)
imgs.volInfo.cluster(1:5, :)

%
%% An example looking at the number of voxels in each voxel list
% Let's take a look at our example:
% Here there are 81,592 voxels in the entire image (E), including "out of mask" voxels
% There are 35,676 "in-mask" voxels (M)
% We start off with 35,676 voxels in the "data list"

% 81592 in the entire image, including "out of mask" voxels
imgs.volInfo.nvox 

% 35676 in the brain mask (with possible data, but some voxels may be empty)
imgs.volInfo.n_inmask

% 35676 voxels in the data list in .dat
size(imgs.dat, 1)  

% 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)  

%
%% 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'));

figure; axis off; montage(mean(imgs), 'trans');

% 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

length(imgs.removed_voxels) % this is still the length of the original mask

% 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
size(imgs.dat, 1) + sum(imgs.removed_voxels) % current voxels + removed must = total number of voxels in mask

%
%% 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');
imgs = resample_space(imgs, atl);

imgs.volInfo.nvox
size(imgs.dat)

% 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}

% 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

% 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.
%