Or, suppose you want to make a parcel-average timecourse using afni's 3dROIstats function, but of a group of parcels together, not each parcel individually. If you make a new mask image with the same number to all the voxels in all the parcels you want treated as a group, 3dROIstats will average them together into a single timecourse.
The final example assigns a new value to each parcel, such as to show the t-values resulting from a statistical test run on each parcel individually.
The code below uses the same files as my NIfTI knitr tutorial, Schaefer2018_400x7_81x96x81.nii.gz and S1200_AverageT1w_81x96x81.nii.gz. To run the code, save the two files somewhere locally, and then change path to their location.
new NIfTI with one parcel only
This first block of code sets the variables for all examples, and then makes a new NIfTI file, p24.nii.gz, which has all the voxels in parcel 24 set to 1, and all the other voxels set to 0.
library(RNifti); # package for NIfTI image reading https://github.com/jonclayden/RNifti
rm(list=ls()); # clear memory
path <- "//storage1.ris.wustl.edu/tbraver/Active/MURI/BraverLab/JoWorking/tutorial/"; # path to input images
p.img <- readNifti(paste0(path, "Schaefer2018_400x7_81x96x81.nii.gz")); # read the parcellation image
# make a new nifti with only parcel 24
new.img <- array(0, dim(p.img)); # blank "brain" of zeros, same size as p.img
new.img[which(p.img == 24)] <- 1; # voxels with value 24 in p.img set to 1
writeNifti(new.img, paste0(path, "p24.nii.gz"), template=paste0(path, "Schaefer2018_400x7_81x96x81.nii.gz"));
# note the template= in the writeNifti function: it specifies that the new file should have the same
# header settings as the original parcellation image; important for both to match properly.
Below are two MRIcron windows, both with S1200_AverageT1w_81x96x81.nii.gz as the anatomical underlay. At left is p24.nii.gz; the title bar text giving the value of the voxel under the crosshairs as 1. The right pane overlay is Schaefer2018_400x7_81x96x81.nii.gz, and the voxel under the crosshairs has the value 24, as expected.
new NIfTI with several parcels given the same value
This code is very similar to the first, but instead of only setting the voxels in parcel 24 to 1, it changes the values of the voxels in all the parcels listed in p.ids to 1 (effectively, making a new parcel out of the previous four):
# make a new nifti with parcels 24, 26, 370, and 231 all given the value 1.
p.ids <- c(24, 26, 231, 360); # parcel ids we want to combine
new.img <- array(0, dim(p.img)); # blank "brain" of zeros, same size as p.img
for (i in 1:length(p.ids)) { new.img[which(p.img == p.ids[i])] <- 1; } # set voxels in p.ids to 1
writeNifti(new.img, paste0(path, "four.nii.gz"), template=paste0(path, "Schaefer2018_400x7_81x96x81.nii.gz"));
new NIfTI with unique numbers for each parcel
To give different numbers to each parcel, loop over all the parcels and assign the corresponding value:
# assign a value to each parcel, such as from the results of a statistical test performed on each parcel individually.
# note: 400 is the total number of parcels in the example parcellation. It is usually best not to hard-code
# the number of parcels, but it is here to keep the code as short as possible.
stat.vals <- rnorm(400); # random numbers to plot, one for each parcel
new.img <- array(0, dim(p.img)); # blank "brain" of zeros, same size as p.img
for (i in 1:400) { new.img[which(p.img == i)] <- stat.vals[i]; } # set voxels in parcel i to stat.value i
writeNifti(new.img, paste0(path, "stat.nii.gz"), template=paste0(path, "Schaefer2018_400x7_81x96x81.nii.gz"));
Below are views of four.nii.gz (left) and stat.nii.gz (right). Since the vector of "statistics" for each parcel in the example code is random, it will be different each time, but the extreme positive and negative values were plotted below in hot and cool colors as is usual for fMRI statistical images.
