I have mutually exclusive binarized masks (around 80) and I want to combine them into one single image where each ROI is represented by a seperate number value in the resulting image - i.e, I want to create an atlas/parcellation with my ROIs.
How do I do this using FSL/AFNI/SPM or other simple commands?
# Let's assume that the images paths are stored in the imgs
from nilearn.image import math_img
new_atlas = math_img('img * 0', img=imgs[0]) # initialize new atlas as image full of 0's
for i, img in enumerate(imgs):
new_atlas = math_img('i1 + (%d + 1) * i2' %i , i1=new_atlas, i2=img)
# With this, value (i+1) in new_atlas corresponds to the ones of the ith image.
In AFNI, there are a variety of ways to make this kind of composite across ROIs.
3dMean. If your data has different values for each ROI - and all the ROIs have no overlap with each other, then you can just use something like
3dMean -max -prefix composite.nii.gz ROI*.nii.gz
where the last part expands to the list of input files. You could also compute the mean, sum or absolute maximum.
3dTstat. This program takes a single dataset as input usually, but you can “trick” it by putting all the files either inside quotes or by making a temporary dataset with all the input datasets. If the ROIs don’t have any particular value (e.g. 1), then you can use the position of the value in a concatenated group of volumes like this:
# make a temporary zero volume with the same grid as the rest
3dcalc -a ROI01+orig. -expr 0 -prefix zero.nii.gz
# chain the zero and the ROIs all together into one dataset
3dTcat -prefix catROIzero.nii.gz zero.nii.gz ROI*.nii.
# find the maximum value position
3dTstat -prefix positionROImax.nii.gz -argmax catROIzero.nii.gz
That puts the position of the maximum in the output dataset, so if the first ROI had a value at a voxel, and all the others did not or if the values were less, then the “argmax” value would be 1 (with the zero placeholder volume). Besides argmax, you could also use max, absmax, median, nzmedian, nzmean, percentile. We use both max and argmax methods to create maximum probability maps from probabilistic atlases where every input volume is a probability or overlap fraction instead of a strict ROI mask. With that fraction, we can also apply a minimum threshold for a combination.
3dcalc. With a list of ROIs, you can also create unique bitwise combinations that account for overlap among ROIs. Here is an example modified a bit from 3dcalc’s help.
0 = none active 1 = A only active 2 = B only active
3 = A and B only 4 = C only active 5 = A and C only
6 = B and C only 7 = all A, B, and C active
3dcalc -a ROI01.nii.gz -b ROI02.nii.gz -c ROI03.nii.gz \
-expr 'step(a)+2*step(b)+4*step(c)' \
-prefix mask_8
This won’t work for 80 regions, but it will be useful for a small number of ROIs.