Creating ggseg-atlas from a ggseg3d-atlas

The best way to create a ggseg-atlas with little manual intervention, is to first make a ggseg3d-atlas, and then convert that to ggseg-atlas. The entire process can take anything from 10minutes, depending on the number of segmentations in the atlas. More segmentations means longer time.

Once you have a ggseg3d-atlas, the function ggseg3d_2_ggseg() can convert that into polygons ready for ggseg.

Prerequisites

dt <- make_ggseg3d_2_ggseg(
  ggseg3d_atlas = dk_3d,
  steps = 1:7, 
  output_dir = tempdir(),
  smoothness = 5,
  tolerance = 0,
  cleanup = FALSE,
  verbose = TRUE)

The function takes few arguments, but does quite a lot in the background. It needs an object that is a ggseg3d-atlas, ad will work with no other input changed. We recommend changing the output_dir to a folder you can keep track of to see how the intermediate results look like. This also makes it possible for you to use the steps argument. If early steps have been run successfully, you can skip them and only run later steps. Re-running steps 5:7 might be necessary to find a good balance between smoothness of the polygons for nicer looking parcellations, and vertex reductions through tolerance to reduce file size. We recommend trying several combinations until you have something that looks nice. Vertex reduction is particularly important, as it significantly reduces the atlas data size, thus making it much faster to plot. This is important for incorporating into Shiny apps, or if you have very many facets to plot. We’d generally recommend a 2d-atlas not have much more than 10k vertices.

Once the data is returned, we recommend doing some data cleaning, particularly of the region column. This column should ideally have human readable names, be short, and have space rather than delimiters. This is to enable nicer labelling when plotting. Once this is done, try plotting it with ggseg.

ggplot() +
  geom_brain(atlas = dt)

DKT example

If you have successfully created a 3d atlas following these steps, you can proceed using the 3d-atlas to make the 2d-atlas. As long as the 3d atlas is correctly made, you should have little difficulty making a 2d version. This pipeline also makes a lot of intermediary files, so be aware of that.

# make atlas ----
dkt <- make_ggseg3d_2_ggseg(dkt_3d,
                            steps = 1:7,
                            smoothness = 2,
                            tolerance = .5,
                            output_dir = "~/Desktop")

# check that it looks ok.
plot(dkt)

Some times, annot-files have a unique colour per hemisphere*region combination. This is does not work well with ggseg. If you fint that the parcellations themselves look fine, but that the colours are off, the atlas palette is surely wrong. You will here need to do some manual editing of the named character vector in the atlas to fix it.

Skipping ggseg3d-atlas

This pipeline is contributed by bbuchsbaum. It is still quite error prone, but was instrumental to creating the above functionality.

The Harvard-Oxford atlas is available with FSL. Here is a process for converting the cortical parts, via FreeSurfer, to a ggseg format.

Converting to FreeSurfer and fixing

  1. Take a copy of the fsaverage folder, and set up SUBJECTS_DIR to point to the folder containing the fsaverage copy.

  2. Import the atlas as a surface

cd ${SUBJECTS_DIR}/fsaverage/mri
mri_vol2surf --sd $(pwd)/../.. --src ${FSLDIR}/data/atlases/HarvardOxford/HarvardOxford-cort-maxprob-thr25-1mm.nii.gz --mni152reg --out  harvard_oxford_surf_rh.mgh --hemi rh --projfrac 0.5

mri_vol2surf --sd $(pwd)/../.. --src ${FSLDIR}/data/atlases/HarvardOxford/HarvardOxford-cort-maxprob-thr25-1mm.nii.gz --mni152reg --out harvard_oxford_surf_lh.mgh --hemi lh --projfrac 0.5
  1. Create a color table. Run an R session in ${SUBJECTS_DIR}/fsaverage/mri
source("ho_ctab.R")
  1. Create a series of individual label files
mkdir Labels
seq 1 48 | parallel mri_vol2label --c harvard_oxford_surf_rh.mgh --id {} --surf fsaverage rh --l ./Labels/rh_HO_{}.label
seq 1 48 | parallel mri_vol2label --c harvard_oxford_surf_lh.mgh --id {} --surf fsaverage lh --l ./Labels/lh_HO_{}.label
  1. Create annotation files
LABS=$(seq 1 48 | parallel echo -n '\ --l ./Labels/lh_HO_{}.label\ ')
mris_label2annot --sd ${SUBJECTS_DIR} --s fsaverage --ctab ../label/ho.annot.ctab ${LABS} --h lh --a ho

LABS=$(seq 1 48 | parallel echo -n '\ --l ./Labels/rh_HO_{}.label\ ')

mris_label2annot --sd ${SUBJECTS_DIR} --s fsaverage --ctab ../label/ho.annot.ctab ${LABS} --h rh --a ho
  1. Check how the surface looks - note lots of holes etc. Using the 0 prob threshold leads to more bleed to neighboring sulci, which I think is harder to fix.
vglrun freeview  --surface ../surf/rh.inflated:annot=../label/rh.ho.annot --surface ../surf/lh.inflated:annot=../label/lh.ho.annot
  1. Fill and smooth
# convert to gifti
mris_convert --annot ../label/rh.ho.annot ../surf/rh.inflated ./fsaverage_rh.label.gii
mris_convert --annot ../label/lh.ho.annot ../surf/lh.inflated ./fsaverage_lh.label.gii

mris_convert --annot ../label/rh.aparc.annot ../surf/rh.inflated ./fsaverage_aparc_rh.label.gii
mris_convert --annot ../label/lh.aparc.annot ../surf/lh.inflated ./fsaverage_aparc_lh.label.gii

./smooth_labels.sh fsaverage_lh_10.label.gii inflated_lh.surf.gii fsaverage_lh_10.smooth.label.gii
./smooth_labels.sh fsaverage_rh_10.label.gii inflated_rh.surf.gii fsaverage_rh_10.smooth.label.gii

## Back to viewing with freesurfer
for hemi in lh rh ; do
mris_convert --annot fsaverage_${hemi}_10.smooth.label.gii inflated_${hemi}.surf.gii ./${hemi}.ho.smooth.annot
done
  1. Screengrabs with tksurfer. This can be reun somewhere other than fsaverage.
./mkPics.sh
  1. Rename the unknowns to medial wall
mv PicsHarvardOxford/lh_\?\?\?_med.tif PicsHarvardOxford/lh_medialwall_med.tif
mv PicsHarvardOxford/rh_\?\?\?_med.tif PicsHarvardOxford/rh_medialwall_med.tif

rm PicsHarvardOxford/*_\?*
  1. Finally, let the R spatial tools work their magic…
source("mkHO.R")

This creates ho_atlas.Rda which contains a couple of data frames, almost ready for inclusion in ggsegExtra.

Comments

I suspect something similar will be possible with tkedit, or other viewers that are able to look at the subcortical structures. Otherwise maximal projections after thresholding will probably suffice.