Combining Atlases Using FSLMaths

Hi Everyone,

Background

I’m using multiple ROIs from a parcelated atlas (Schaefer et al., 2018) to run an analysis. It has 400 parcels but I’m hoping to generate a new atlas excluding any parcels or ROIs that are not in the Default Mode Network (DMN). I don’t want a single mask covering all of the DMN regions, but rather an atlas that clearly delineates between each of the DMN ROIs. I’m using thresholding via fslmaths to eliminate the ROIs that are not part of the DMN. For example:

# Atlas path 
ATLAS=/data/tools/schaefer_parcellations/MNI/
SCRIPTS=/data/project/scripts

# Using fslmaths to exclude any ROI whose index or value is greater than 41
# ROIs 1 through 41 contain DMN-related regions in the left hemisphere
fslmaths ${ATLAS}/Schaefer2018_400Parcels_Kong2022_17Networks_order_FSLMNI152_2mm.nii.gz \
         -uthr 41 \
         ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_L_DMN.nii.gz

# DMN-related regions in the right hemisphere correspond to ROI indices 201 through 236.
# As such I'll use thresholding again to generate an atlas containing only those regions 
fslmaths ${ATLAS}/Schaefer2018_400Parcels_Kong2022_17Networks_order_FSLMNI152_2mm.nii.gz \
         -thr 201 \
         ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz

fslmaths ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz \
         -uthr 236 \
         ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz

This leaves me with a separate atlas for left and right cortices which I need to recombine. Before I do so, I want to reindex the right hemisphere ROIs to be 42 through 77 rather than 201 through 236. This way, when I combine them the ROIs indices will be 1 through 77. I use the subtract argument for that:

# Reindexing values from the right hemisphere
fslmaths ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz \
         -sub 159 \
         ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz

This appears to work great. When I view the two separate atlases the left ROI indices range from 1 to 41 and the right indices range from 42 to 77. I then combine the two using the add function:

# Combining the Atlases
fslmaths ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_L_DMN.nii.gz \
         -add ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz \
         ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_DMN.nii.gz

The Problem:

However, after doing this I’m left with this issue wherein the ROIs in the right cortex remain indexed as 42 through 77 and the ROIs on the left hemisphere are re-indexed as -159 through -118. It seems like after being recombined they somehow redefined their position according to how far away they were from the greatest index. (i.e., Because the highest value index changed from 236 to 77, a change of 159, all left hemisphere ROIs also were reduced by 159).

Am I misunderstanding what the -add, -thr, and -sub arguments are fundamentally doing? Is there perhaps an easier way to do what I’m trying to do? Of course, if any of the code or details above are not reproducing well or are confusing, please let me know. I appreciate the help!

Hi,

One issue is coming from this line:

# Reindexing values from the right hemisphere
fslmaths ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz \
         -sub 159 \
         ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz

Indeed, it would change your labels ranging from 201 to 236 to labels ranging from 42 to 77. However, the rest of the voxels, whose intensity was set to zero after thresholding, is now at -159!

You just need to add a new line after the command above:

 fslmaths ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz \ 
               -thr 0 \
               ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz

and go on with the rest of your script.

Also:

# DMN-related regions in the right hemisphere correspond to ROI indices 201 through 236.
# As such I'll use thresholding again to generate an atlas containing only those regions 
fslmaths ${ATLAS}/Schaefer2018_400Parcels_Kong2022_17Networks_order_FSLMNI152_2mm.nii.gz \
         -thr 201 \
         ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz

fslmaths ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz \
         -uthr 236 \
         ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz

Could just be written in one command:

fslmaths ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz \
         -uthr 236 -thr 201 \
         ${SCRIPTS}/dir_ROIs/Schaefer_400Parcel_R_DMN.nii.gz
1 Like

Oh man, of course! Seems so obvious in retrospect. I spent a couple hours scratching my head over this and reapplying -thr did not occur to me. Thank you so much. I really appreciate it.

Hiya! This thread is exactly what I need - if I only the solution by @jsein worked as well for me as it did for @wjpmitchell3! Would either of you mind taking a look to see if I’ve implemented it incorrectly?

I have two atlases which I would like to merge. The first is the 200 region Schaefer atlas, and the second, and the Scale 1 version of the Melbourne subcortical atlas (paritioned into 16 regions). They’ve both been registered to standard 2mm MNI space.

Since the Schaefer atlas has 200 regions, I figured I’d add 200 to the Melbourne atlas, so the final atlas ranged from 1-216. However, when I employ the following code, the intensity values range from 0-200. Oddly, when I remove the middle line, I get the same result. Any advice is welcome, and thank you for your time!

fslmaths Melbourne.nii -add 200 Melbourne.nii
fslmaths Melbourne.nii -thr 0 Melbourne.nii
fslmaths Melbourne.nii -add Schaefer.nii Combo.nii.gz

Hi @daniellekurtin, welcome to Neurostars!

I am surprised that it doesn’t work, I did a try on my side and the same commands seem to achieve the desired result.
However, as you noticed the second line is useless here as all the values from Melbourne.nii are already above 0 (and above 200 at this stage to be more precise).

Did you check that Melbourne.nii has indeed values between 200 and 216 after the first fslmaths command?
Are you sure that Combo.nii.gz has not value above 200?
As a side note, it may be more convenient to give a different name to the resulting image after fslmaths, especially during the debugging stage, to not overwrite the original file.

Heyya @jsein, thanks for the welcome! :slight_smile:
Yep, using “fslstats Melbourne.nii -R” returns min-max intensities of 0-16 before adding 200, and after, the min-max is 200-216 (which made me realise - I should be adding 201 to Melbourne, rather than 200!)
Ahaha good pickup on the same name after using fslmaths; originally I had used different names, but since @wjpmitchell3 used the same name after -thr, -sub, and -uthr functions, I kept with the thread’s convention.
I’m not sure what else to do. Since I’m trying to combine the atlases so it streamlines an fslmeants timeseries extraction, my worst-case patch is the extract timeseries from each atlas, then combine the outputs. And that’s not too bad of a worse-case scenario.

I looked a bit more and found out two things:

  • There was one line missing (it works if you keep the same filename all around but as I said, for testing it may be more convenient to use intermediate names):
fslmaths Melbourne.nii -add 200 Melbourne200.nii
fslmaths Melbourne200.nii -thr 201 Melbourne200thr.nii
fslmaths Schaefer.nii -add Melbourne200thr.nii Combo.nii.gz 

=> That way you avoid to have a background with values at 200!

  • With fslstats, the result is
fslstats Combo.nii.gz -R
0.000000 396.000000

and you can see the result on the capture below:

You can see the value of the voxel when the green crosshair is at for the different images in the right bottom corner.

Also not that in my case there was some overlap between the two atlas, that is why the maximum of the combined atlas is at 396.

Tell us if you are not able to reproduce this.

Also, another thing: In my case, even if Melbourne.nii and Schaefer.nii are in the same MNI space, on is store in LAS format and the other is stored in RAS. FSL and FSLeyes are clever about it but it is important to check that out and preferably put them in the same convention.

Heyya @jsein, thank you SO MUCH for your suggested solution - it absolutely worked! I was not aware of the RAS and LAS formats for storage, so that’s my next reading rabbit hole to dive down - thank you for bringing that to my attention :slight_smile: