How to Get Grey Matter Volume of Specific Brain Regions in SPM/MATLAB?

I want to extract grey matter volume for specific regions (e.g., prefrontal cortex, amygdala) from MRI data that has been preprocessed. My lab gave me a code (see below) for calculating global grey matter volume (and white matter volume and CSF), but I am not sure whether I can adapt this script to also extract regional grey matter volume. Would appreciate help!

% Get total volumes in mm3

grey_totals(i) = get_totals ([data_path num2str(subnum) '/c1' num2str(subnum) '.nii']);

white_totals(i) = get_totals ([data_path num2str(subnum) '/c2' num2str(subnum) '.nii']);

csf_totals(i) = get_totals ([data_path num2str(subnum) '/c3' num2str(subnum) '.nii']);

end

[Karl]

It would be useful to have some additional information:

1. How is the get_totals() function defined? In particular, what does it expect as input?
2. In the example, does each .nii file correspond to a binary mask representing the named material?
3. For each specific region in which you're interested, do you have an .nii file corresponding to a binary mask representing the region?

[Arjun]

Of course.

1. function [t, files] = get_totals(files, thr, msk)
2. In the example, the .nii files correspond to the segmented tissues (i.e., grey matter, white matter, and CSF) for each subject
3. I do have masks for each ROI that I created from wfu_pickatlas

Thank you so much for your help!

[Karl]

Thanks for the additional information. I think that you should be able to adapt your example script to extract the regional greyscale matter volumes. I'd suggest trying first for a single region for a single subject, something like:

region_totals = get_totals([path_to_ROI_mask]);

where you set path_to_ROI_mask to be the path to one of the masks that you've created with WFU_PickAtlas.

If the above doesn't work, or gives an answer that's obviously wrong, could you upload an example mask file? If it does work, it should then just be a case of naming your mask files, and organising in folders, so that you can process masks for all subjects in a loop. (I'm guessing that i in the example is a loop index.)

[Arjun]

I ran the code below:

>> %get SPM12
addpath '[path to SPM12]';
spm
>> get_region = get_totals(['[path to ROI mask']);
>> get_region(['[path to MRI data]'  num2str(subnum) '/c1' num2str(subnum) '.nii']);

But I got the following error message:

Index exceeds the number of array elements. Index must not exceed 1.

I am new to coding, so I am not too sure what this means. I am also uploading the mask file that I tried using, if that helps.

Thank you so much again for your help!

[Karl]

Thanks for uploading the example NIfTI mask. After unzipping, the code:

info = niftiinfo('Amygdala_mask.nii');
mask = niftiread(info);
volshow(mask);

gives:

and the volume in cubic millimetres is:

volume = sum(mask(:)) * prod(info.PixelDimensions)

This is one way of obtaining the region volume, but using get_totals() should also work. In your code, the path to the data doesn't look right. If the (unzipped) mask file were in the folder from which you run MATLAB, you should be able to use:

amygdala_totals = get_totals('Amygdala_mask.nii');

In general, you might not be running MATLAB from the folder where you have your data, so you'd need something like:

amygdala_totals = get_totals('/path/to/my/data/Amygdala_mask.nii');

In case useful, there's some information about file paths in the documentation:

https://mathworks.com/help/matlab/matlab_env/specify-file-names.html

[Arjun]

In general, you might not be running MATLAB from the folder where you have your data, so you'd need something like:

amygdala_totals = get_totals('/path/to/my/data/Amygdala_mask.nii');

With reference to the comment above, I am currently running the code:

%% get amygdala grey, white matter, and CSF volumes

%get SPM12
addpath '/Users/arjunthanaraju/Downloads/PhD/spm12/';
spm

for i = 1:length(files_mri)

subnum = files_mri(i).name; %get names of each subject's folder

%add amygdala mask to get_totals function
amygdala_totals = get_totals(['/Volumes/SEAGATE_BAC/Amygdala_Mask.nii']);

%get amygdala volumes in mm3
amygdala_grey_totals(i) = amygdala_totals (['/Volumes/SEAGATE_BAC/397Participants/'  num2str(subnum) '/c1' num2str(subnum) '.nii']);
amygdala_white_totals(i) = amygdala_totals (['/Volumes/SEAGATE_BAC/397Participants/'  num2str(subnum) '/c2' num2str(subnum) '.nii']);
amygdala_csf_totals(i) = amygdala_totals (['/Volumes/SEAGATE_BAC/397Participants/'  num2str(subnum) '/c3' num2str(subnum) '.nii']);
end 

%combine into table
globals = table(subjects.subjects, amygdala_grey_totals', amygdala_white_totals', amygdala_csf_totals');

%convert Var1 column name to Subjects
globals.Properties.VariableNames{'Var1'} = 'Subjects';

When I run this code, I get the previously stated error:

Index exceeds the number of array elements. Index must not exceed 1.

Can you point me towards where I made a mistake in the code above? Thank you so much!

[Karl]

The error comes from trying to access amygdala_totals as anything other than a scalar - the same as the following:

amygdala_totals = 3.744;
x = amygdala_totals('path');

Index exceeds the number of array elements. Index must not exceed 1.

Could you try running just the following:

%get SPM12
addpath '/Users/arjunthanaraju/Downloads/PhD/spm12/';

%obtain volume of amygdala
amygdala_totals = get_totals(['/Volumes/SEAGATE_BAC/Amygdala_Mask.nii']);
disp(['Volume of amygdala = ', num2str(amygdala_totals), ' (units to be checked)']);

This should hopefully confirm that the basic volume calculation is working for you.

[Arjun]

Could you try running just the following:

%get SPM12
addpath '/Users/arjunthanaraju/Downloads/PhD/spm12/';
%obtain volume of amygdala
amygdala_totals = get_totals(['/Volumes/SEAGATE_BAC/Amygdala_Mask.nii']);
disp(['Volume of amygdala = ', num2str(amygdala_totals), ' (units to be checked)']);

I could run this code without issue! What would be next?

[Karl]

Excellent! A possible next step would be to organise the ROI masks for which you want to extract volume information in a way that makes them easy to analyse.

In the example code that you have, there's a folder for each subject, and each subject folder contains multiple ROI masks. You could do something similar. For example, if the first subject has identifier 12345, you could create a folder /Volumes/SEAGATE_BAC/397Participants/12345, and inside that place the ROI masks, for example amygdala_12345.nii, some_other_roi_12345.nii.

You could then loop over subjects, extracting volumes for each ROI. If in your code:

subnum = files_mri(i).name;

correctly defines the subject identifier, you should be able to obtain volumes with lines like:

amygdala_totals(i) = get_totals(['/Volumes/SEAGATE_BAC/397Participants/'  num2str(subnum) '/amygdala_' num2str(subnum) '.nii']);
some_other_roi_totals(i) = get_totals(['/Volumes/SEAGATE_BAC/397Participants/'  num2str(subnum) '/some_other_roi_' num2str(subnum) '.nii']);

You'd then change the definition of the globals table to reflect the ROIs for which you have data.

[Arjun]

Hi! Unfortunately, I am still unable to get the code to work for me but I suspect this is an issue with the looping rather than the actual get_totals function. I will continue working on this but thank you so much for your kind help! Your guidance and explanations have been very useful. I really appreciate it!

[Karl]

You're welcome! If you're unable to get the looping to work, please do submit another question, ideally including code and files for reproducing the problem (as you did in the comments here). Hopefully I, or someone else, will be able to help!

[Arjun]

Thank you so much! I definitely will.

Answer 1

[Karl]

Summarising from the comments, if 'roi.nii' is a NIfTI file containing a binary mask representing a region of interest (ROI), and the Statistical Parametric Mapping (SPM12) package is on the MATLAB path, the ROI volume can be calculated using the get_totals() function written by Ged Ridgway:

% Calculate ROI volume using SPM12.
roi_volume = get_totals('roi.nii');

The ROI volume can also be calculated without having SPM12 installed:

% Write example ROI (cuboid) in NIfTI format.
mask = zeros(100,100,100);
[x1, y1, z1] = deal(31, 41, 11);
[dx, dy, dz] = deal(40, 20, 80);
mask(x1:x1+dx-1,y1:y1+dy-1,z1:z1+dz-1) = 1;
niftiwrite(mask,'roi.nii')

% Read ROI from NIfTI file, and calculate its volume.
roi_info = niftiinfo('roi.nii');
roi = niftiread(roi_info);
roi_volume = sum(roi(:)) * prod(roi_info.PixelDimensions)

roi_volume = 64000

% Check that the calculated volume is as expected from the cuboid dimensions.
assert(roi_volume == dx * dy * dz)

Whichever method is used, the volume units should be checked.

[Arjun]

Just wanted to add on for future users that the CAT12 toolbox in SPM also provides estimates for each ROI based on selected atlases. This made it very easy for me to extract ROI values for both grey matter volume and cortical thickness!

You can read more in the CAT12 manual available here: https://neuro-jena.github.io/cat12-help/