Automatically masking the medial wall when accessing morphometry data

Typically the vertices of the medial wall are excluded from the analysis, as they are not part of the cortex. If you have the files label/lh.cortex.label and label/rh.cortex.label available for your subjects, you can automatically mask non-cortex data. The morhometry values for the vertices outside the cortex will then be set to NA. To do this, pass the parameter cortex_only = TRUE to the morphometry data functions. Here is an example:

# Load the full hemisphere data, including media wall:
fulldata = group.morph.native(subjects_dir, subjects_list, "thickness", "lh");
mean(fulldata$subject1);

# This time, restrict data to the cerebral cortex:
cortexdata = group.morph.native(subjects_dir, subjects_list, "thickness", "lh", cortex_only=TRUE);
mean(cortexdata$subject1, na.rm=TRUE);

This also works for standard space data. Keep in mind that in that case, you need to have the files label/lh.cortex.label and label/rh.cortex.label available for your template subject.

Labels and masks

You can turn a label (or the combination of several labels) into a mask. A mask for a surface is a logical vector that contains one value for each vertex in the surface. Masks are very convenient for selecting subsets of your data for analysis. Here is an example that first creates a label from a region of an annotation, then it creates a mask from the label. Finally, it modifies the existing mask by adding more vertices from a second region:

surface = 'white';
hemi = 'both';
atlas = 'aparc';
region = 'bankssts';

# Create a mask from a region of an annotation:
lh_annot = subject.annot(subjects_dir, subject_id, 'lh', atlas);
rh_annot = subject.annot(subjects_dir, subject_id, 'rh', atlas);
lh_label = label.from.annotdata(lh_annot, region);
rh_label = label.from.annotdata(rh_annot, region);
lh_mask = mask.from.labeldata.for.hemi(lh_label, length(lh_annot$vertices));
rh_mask = mask.from.labeldata.for.hemi(rh_label, length(rh_annot$vertices));

# Edit the mask: add the vertices from another region to it:
region2 = 'medialorbitofrontal';
lh_label2 = label.from.annotdata(lh_annot, region2);
rh_label2 = label.from.annotdata(rh_annot, region2);
lh_mask2 = mask.from.labeldata.for.hemi(lh_label2, length(lh_annot$vertices),
  existing_mask = lh_mask);
rh_mask2 = mask.from.labeldata.for.hemi(rh_label2, length(rh_annot$vertices),
  existing_mask = rh_mask);

A single morphometry measure in native space

The fsbrain package provides a high-level interface to access data for groups of subjects. Here, we compute the mean thickness for each subject in native space:

mean_thickness_lh_native = group.morph.agg.native(subjects_dir, subjects_list, "thickness", "lh", agg_fun=mean);
mean_thickness_lh_native;
# output:
#  subject_id hemi measure_name measure_value
#1   subject1   lh    thickness      2.437466
#2   subject2   lh    thickness      2.437466

Note that subject1 and subject2 in the test data are identical (2 is a copy of 1), so we get the same results for both.

A single morphometry measure in standard space

Here, we compute the mean thickness for each subject in standard space:

mean_thickness_lh_std = group.morph.agg.standard(subjects_dir, subjects_list, "thickness", "lh", fwhm="10", agg_fun=mean);
mean_thickness_lh_std;
# output:
  subject_id hemi measure_name measure_value
1   subject1   lh    thickness       2.32443
2   subject2   lh    thickness       2.32443

Several measures over several hemispheres in native space

When working with native space data for groups, vertex-wise comparison is not possible and one is often interested in aggregating data or summary statistics. Obtaining these for a group becomes easy with ‘fsbrain’:

agg_nat = group.multimorph.agg.native(subjects_dir, subjects_list, c("thickness", "area"), c("lh", "rh"), agg_fun = mean);
head(agg_nat);
# output:
#  subject_id hemi measure_name measure_value
#1   subject1   lh    thickness     2.4374657
#2   subject2   lh    thickness     2.4374657
#3   subject1   lh         area     0.6690556
#4   subject2   lh         area     0.6690556
#5   subject1   rh    thickness     2.4143047
#6   subject2   rh    thickness     2.4143047

The measure values are the mean thickness/area/volume over all vertices of the subject, computed from the respective native space morphometry files (e.g., ‘subject1/surf/lh.thickness’, ‘subject2/surf/lh.thickness’, and so on).

Note that you can get short format by setting ‘cast = FALSE’, which will result in a dataframe with the following columns instead:

agg_nat2 = group.multimorph.agg.native(subjects_dir, subjects_list, c("thickness", "area"), c("lh", "rh"), agg_fun = mean, cast=FALSE);
head(agg_nat2);
# output:
# subject_id lh.thickness   lh.area rh.thickness   rh.area
#1   subject1     2.437466 0.6690556     2.414305 0.6607554
#2   subject2     2.437466 0.6690556     2.414305 0.6607554

Several measures over several hemispheres in standard space

You could do the same in standard space with group.multimorph.agg.standard(), even though it’s less common as typically vertex-wise analysis is used in standard space. This requires that you pass the fwhm parameter to identify which smoothing value you want:

data_std = group.multimorph.agg.standard(subjects_dir, subjects_list, c("thickness", "area"), c("lh", "rh"), fwhm='10', agg_fun = mean);
head(data_std);
# output:
#  subject_id hemi measure_name measure_value
#1   subject1   lh    thickness     2.3244303
#2   subject2   lh    thickness     2.3244303
#3   subject1   lh         area     0.5699257
#4   subject2   lh         area     0.5699257
#5   subject1   rh    thickness     2.2926377
#6   subject2   rh    thickness     2.2926377

Other parameters exist to define the template subject, which defaults to ‘fsaverage’ in the example above.

Automatically masking the medial wall when aggregating

For the hemisphere level, one can also automatically mask data from the medial wall, as explained above. The only thing to keep in mind is that your aggregation function must be able to cope with NA values in this case. If the function requires a special parameter to be passed to it for this to work out, use the agg_fun_extra_params parameter. Here is an example with mean as the aggregation function:

data_std_cortex = group.multimorph.agg.standard(subjects_dir, subjects_list, c("thickness", "area"), c("lh", "rh"), fwhm='10', agg_fun = mean, cortex_only=TRUE,  agg_fun_extra_params=list("na.rm"=TRUE));
head(data_std_cortex);

This will ensure that mean gets passed the extra parameter na.rm=TRUE. You will have to adapt this based on the aggregation function.

Native space

This is typically done in native space, and here is an example:

atlas = 'aparc';         # or 'aparc.a2009s', or 'aparc.DKTatlas'.
measure = 'thickness';
region_means_native = group.agg.atlas.native(subjects_dir, subjects_list, measure, "lh", atlas, agg_fun = mean);
head(region_means_native[,1:6]);
# output:
#          subject bankssts caudalanteriorcingulate caudalmiddlefrontal   cuneus entorhinal
#subject1 subject1 2.485596                 2.70373            2.591197 1.986978   3.702725
#subject2 subject2 2.485596                 2.70373            2.591197 1.986978   3.702725

We only showed the first 6 columns here.

The measure values are the mean thickness over all vertices of the respective region of the subject, computed from the respective native space morphometry files (e.g., ‘subject1/surf/lh.thickness’, ‘subject2/surf/lh.thickness’, and so on). The brain parcellation for the subject is read from its annotation file ( ‘subject1/label/lh.aparc.annot’ in this case).

Standard space

You could do the same in standard space with group.agg.atlas.standard(), even though it’s less common as typically vertex-wise analysis is used in standard space. This requires that you pass the fwhm parameter to identify which smoothing value you want, and that you have the template subject in the subjects_dir. The template subject defaults to fsaverage.

region_means_std = group.agg.atlas.standard(subjects_dir, subjects_list, measure, "lh", atlas, fwhm = '10', agg_fun = mean);
head(region_means_std[1:5]);
# output:
#          subject bankssts caudalanteriorcingulate caudalmiddlefrontal   cuneus
#subject1 subject1 2.583408                2.780666            2.594696 2.018783
#subject2 subject2 2.583408                2.780666            2.594696 2.018783

Other parameters exist to define the template subject, which defaults to ‘fsaverage’ in the example above.

Visualizing annotations (brain parcellations based on an atlas)

Let’s first visualize an annotation:

vis.subject.annot(subjects_dir, 'subject1', 'aparc', 'both', views=c('si'));

This will give you an interactive window in which you can freely rotate a view similar to the one shown in the following screenshot:

Figure 1: Atlas regions visualized on a brain surface mesh.

Visualizing morphometry data

You can also visualize morphometry data:

vis.subject.morph.native(subjects_dir, 'subject', 'thickness', hemi='both', views=c('si'))

You want want to load the data first if you want to do advanced pre-processing or other computations before visualizing it:

morph_data_lh = subject.morph.native(subjects_dir, 'subject1', 'thickness', 'lh');
morph_data_rh = subject.morph.native(subjects_dir, 'subject1', 'thickness', 'rh');
# Do something with the morph_data here.
vis.data.on.subject(subjects_dir, 'subject1', morph_data_lh, morph_data_rh, views=c('si'));

Figure 2: Cortical thickness data visualized on a brain mesh in interactive view.

Visualizing labels

To visualize a label, use the vis.subject.label function:

surface = 'white';
hemi = 'both';
label = 'cortex.label';
vis.subject.label(subjects_dir, subject_id, label, hemi);

surface = 'white';  # If possible, use the 'inflated' surface instead: it is much easier to find the vertices on it. We do not
#  use it here because the inflated surface is not shipped with the example data for this package to reduce download size.

# For the left hemi, we just specify 3 vertices. They are very small in the high-resolution mesh and may be hard to spot.
my_lh_result_cluster_vertices = c(1000, 1001, 1002);  # Would be a lot more than just 3 in real life, and they would be adjacent (these are not).
lh_labeldata = my_lh_result_cluster_vertices;

# For the right hemi, we extend the neighborhood around our vertices of interest for the visualization. This makes the a lot easier to spot.
my_highest_effect_size_vertex = 5000;
rh_labeldata = c(my_highest_effect_size_vertex);
rh_surface = subject.surface(subjects_dir, subject_id, surface, 'rh');
rh_labeldata_neighborhood = mesh.vertex.neighbors(rh_surface, rh_labeldata);   # extend neighborhood for better visibility
rh_labeldata_neighborhood = mesh.vertex.neighbors(rh_surface, rh_labeldata_neighborhood$vertices);   # extend neighborhood again

# Hint: Check the area around the visual cortex when searching for the vertices in interactive mode.
vis.labeldata.on.subject(subjects_dir, subject_id, lh_labeldata, rh_labeldata_neighborhood$vertices, views=c('si'), surface=surface);

Visualizing masks

In the following example, we create a mask by combining several labels and visualize it:

surface = 'white';
hemi = 'both';
atlas = 'aparc';
region = 'bankssts';

# Create a mask from a region of an annotation:
lh_annot = subject.annot(subjects_dir, subject_id, 'lh', atlas);
rh_annot = subject.annot(subjects_dir, subject_id, 'rh', atlas);
lh_label = label.from.annotdata(lh_annot, region);
rh_label = label.from.annotdata(rh_annot, region);
lh_mask = mask.from.labeldata.for.hemi(lh_label, length(lh_annot$vertices));
rh_mask = mask.from.labeldata.for.hemi(rh_label, length(rh_annot$vertices));

# visualize it
vis.mask.on.subject(subjects_dir, subject_id, lh_mask, rh_mask);

Visualizing results of region-based analyses

The result of a region-based analysis often is one result value (e.g., a p-value or an effect size) per atlas region. In the following example, we visualize such data on the fsaverage template subject, but you can use any subject of course.

atlas = 'aparc';
template_subject = 'fsaverage';
# Some directory where we can find the template_subject. This can be omitted if FREESURFER_HOME or SUBJECTS_DIR is set and the template subject is in one of them. In that case, the function will find fsaverage in there by default. Also see the function download_fsaverage().
template_subjects_dir = "~/software/freesurfer/subjects";    # adapt to your machine


# For the left hemi, we manually set data values for some regions.
lh_region_value_list = list("bankssts"=0.9, "precuneus"=0.7, "postcentral"=0.8, "lingual"=0.6);

# For the right hemisphere, we do something a little bit more complex: first get all atlas region names:
atlas_region_names = get.atlas.region.names(atlas, template_subjects_dir=template_subjects_dir, template_subject=template_subject);
# As mentioned above, if you have fsaverage in your SUBJECTS_DIR or FREESURFER_HOME is set, you could replace the last line with:
#atlas_region_names = get.atlas.region.names(atlas);

# OK, now that we can check all region names. We will now assign a random value to each region:
rh_region_value_list = rnorm(length(atlas_region_names), 3.0, 1.0);         # create 36 random values with mean 3 and stddev 1
names(rh_region_value_list) = atlas_region_names;                           # use the region names we retrieved earlier

# Now we have region_value_lists for both hemispheres. Time to visualize them:
vis.region.values.on.subject(template_subjects_dir, template_subject, atlas, lh_region_value_list, rh_region_value_list);

See the help for the vis.region.values.on.subject function for more options, including the option to define the default value for the regions which do not appear in the given lists.

Figure 3: Various fsbrain visualizations. A Native space morphometry data. B Results of a region-based analysis. C An atlas.

Full control: visualizing vertex colors and color layers

You can visualize pre-defined vertex colors using the vis.color.on.subject function. All you need is one vector of colors (one color per vertex) for each hemisphere. You can create the vectors yourself, or use the intermediate level API of fsbrain to create such a color layer from various types of data. Here is an example that loads a colorlayer from mean curvature morphometry data, manipulates the color values (desaturates them, i.e., turns them into grayscale), and visualizes the result:

# load mean curv data
meancurv = subject.morph.native(subjects_dir, subject_id, "curv", "both", split_by_hemi = TRUE);

# curvature data is noisy, clip it for better visualization (optional).
meancurv = lapply(meancurv, clip.data);

# desaturate colors for left hemisphere
cl$lh = desaturate(cl$lh);

# visualize it
vis.color.on.subject(subjects_dir, subject_id, cl$lh, cl$rh);

Technically, a color layer is just a nemd list with entries ‘lh’ and/or ‘rh’, each of which is a vector of colors.

Figure size

The default output resolution (device window size) in rgl is very low, and it is highly recommended to increase it. The default size may be so small that there is not enough space to draw a colorbar. You can set the device window size globally. E.g., if you want the window to have 1200 x 1200 pixels, and the window should be placed at screen pixel position (50, 50), i.e., near the upper left corner of your screen, you would do:

library('rgl');
r3dDefaults$windowRect=c(50, 50, 1200, 1200);

Now call any fsbrain visualization function, and the setting will be used.

It is also possible to set this when calling an fsbrain visualization function, in which case it is only applied to the current function call. For any visualization function that accepts a parameter named ‘rgloptions’, use something like this:

rgloptions = list('windowRect'=c(50, 50, 1200, 1200));

Different views

You can visualize the data in different views. So far, we have used only the si view. The following views are available:

• si or single interactive view: renders a single instance of the mesh and allows the user to interactively explore the scene (rotate/zoom with the mouse). This view is good for live inspection of data during analysis.
• sr or single rotating view: renders a single instance of the mesh rotates the camera around the mesh. Can be combined with rglactions (see section ‘Creating movies’ below) to create animated GIF images. This view is good for live inspection of data, and the resulting movies can be used in presentations.
• t4 or tiled, 2x2, 4 angles view: renders the meshes from four different angles. This kind of view is often used in publications.
• t9 or tiled, 3x3, 8 angles view: renders the meshes from eight different angles. This kind of view is often used in publications.

To change the view, set the view parameter in any visualization function. You can specify more than one view if needed:

vis.subject.morph.native(subjects_dir, 'subject', 'thickness', hemi='both', views=c('t4', 't9'))

See Figure 3A for an example of view t9 and Figure 3B for an example of view t4.

Changing the visualisation surface

By default, the data is visualized on the white surface. For some data, other surfaces may be more apppropriate. You can visualize the data on any FreeSurfer surface you have available. To do this, just set the surface parameter in any visualization function. The most commonly used surfaces and white, pial and inflated. Examples:

vis.subject.morph.native(subjects_dir, 'subject', 'thickness', hemi='both', views=c('si'), surface='inflated')

See Figure 4 for examples for the different surfaces.

Figure 4: Different visualization surfaces. A white surface. B pial surface. C inflated surface.

Changing the resolution of figures (or other rgloptions)

The resolution of rgl windows is set in the call to rgl::par3d. You can pass arbitrary options to the call by specifying the parameter rgloptions when calling any visualization function. Here is an example that increases the resolution of the output window to 800x800 pixels and opens the window at screen position 50, 50:

rgloptions = list("windowRect"=c(50, 50, 1000, 1000));
vis.subject.morph.native(subjects_dir, 'subject', 'thickness', hemi='both', views=c('si'), rgloptions=rgloptions)

See the documentation of rgl::par3d for all available options.

Note that for some OpenGL implementations, it may lead to artifacts or other problems if the window is not fully visible while actions like taking a screenshot or creating a movie are performed. In such a case you may have to adapt the position of the window on the screen (the first 2 parameters of the windowRect setting) to make it fully visible if parts of it are off-screen. This also means you should not perform other actions on the machine while a movie is being recorded, as taking away the focus from the rgl window or opening other windows on top of it may lead to artifacts.

Using backgrounds when visualizing data

By default, when you plot vertex-wise data that contains NA values, the respective parts of the surface are rendered in white (cf. Figure 3B), or whatever color your colormap function assigns to NA values. One can set a different background using the bg option of the vis functions. Typical backgrounds are binarized gray-scale visualizations of mean curvature or sulcal depth, as they show the pattern of gyri and sulci and thus make it easier to see where the clusters are, especially with the inflated surface. Here is an example that plots clusters against a background of mean curvature:

subjects_dir = find.subjectsdir.of("fsaverage")$found_at;
subject_id = 'fsaverage';

lh_demo_cluster_file = system.file("extdata", "lh.clusters_fsaverage.mgz", package = "fsbrain", mustWork = TRUE);
rh_demo_cluster_file = system.file("extdata", "rh.clusters_fsaverage.mgz", package = "fsbrain", mustWork = TRUE);
lh_clust = freesurferformats::read.fs.morph(lh_demo_cluster_file);   # contains a single positive cluster
rh_clust = freesurferformats::read.fs.morph(rh_demo_cluster_file);   # contains two negative clusters
vis.symmetric.data.on.subject(subjects_dir, subject_id, lh_clust, rh_clust, bg="sulc");

This example requires the fsaverage subject. Note that the vis.symmetric.data.on.subject function by default maps zero values to NA (and uses a diverging colormap with symmetric borders and a transparent color for NA data values), so you do not have to do anything more in this case. The output is displayed in Figure 5A.

Figure 5: Different backgrounds and color layers. A Clusters on white fsaverage surface with sulc background. B Region-wise p-values with curv background, inflated surface. C Color layer displaying outlines of aparc atlas regions in the respective colors.

The following backgrounds are available:

• ‘curv’: binarized mean curvature, gray-scale. Based on the lh.curv and/or rh.curv files for the subject.
• ‘sulc’: binarized sulcal depths, gray-scale. Based on the lh.sulc and/or rh.sulc files for the subject.
• ‘aparc’: desatured aparc atlas. Based on the lh.aparc.annot and/or rh.aparc.annot files for the subject.

Technically, the backgrounds are color layers, and you are free to create your own and pass it as the value for the bg parameter instead of using one of the pre-defined character strings listed above. Alpha blending is used to merge the different layers into the final colors.

Keep in mind that the background will be visible only through the (semi-)transparent parts of your foreground data. Morphometry or atlas data typically contains a fully opaque color value for each vertex, and thus, the background will not be visible. You can manually add transparency to (parts of) your color data by using color layers as explained above.

Saving screenshots and creating GIF movies (and other rglactions)

The rglactions parameter can be passed to any visualization function to trigger certain actions during the visualization. These are actions like saving a screenshot or a movie or animation in GIF format. Some rglactions only make sense in certain views (e.g., recording a movie only makes sense in animated views, i.e., the sr view). Views will silently ignore actions which do not make sense for them.

The following rglactions are available:

• Take a screenshot in PNG format. Key: “snapshot_png”. Value: string, the path of the output file, including the file extension. Path expansion is performed. Example: rglactions = list("snapshot_png"="~/fsbrain.png").
• Record a movie. Only handled in rotating views like sr. Key: “movie”. Value: string, the base filename of the movie. The movie will be saved in your home directory, with file extension gif. Example: rglactions = list("movie"="fsbrain_rotating_brain")
• Clip the data to remove outliers before visualization. Data smaller than the lower percentile and larger than the upper percentile will be set to the respective percentile value. Required to visualize data with outliers or extreme values, like mean or Gaussian curvature. Key: “clip_data”. Value: vector of length 2 with the lower and upper percentile. Example: rglactions = list("clip_data"=c(0.05, 0.95))

The actions can be combined, e.g., to clip the data and take a screenshot: rglactions = list("clip_data"=c(0.05, 0.95), "snapshot_png"="~/fsbrain.png").

subjects_dir = fsbrain::get_optional_data_filepath("subjects_dir");
subject_id = 'subject1';
rgloptions=list("windowRect"=c(50, 50, 600, 600));     # the first 2 entries give the position on screen, the rest defines resolution as width, height in px
surface = 'white';
measure = 'thickness';
movie_base_filename = sprintf("fsbrain_%s_%s_%s", subject_id, measure, surface);
rglactions = list("movie"=movie_base_filename, "clip_data"=c(0.05, 0.95));
# Creating a movie requires the rotating view ('sr' for 'single rotating'). The action will be silently ignored in all other views.
vis.subject.morph.native(subjects_dir, subject_id, measure, 'both', views=c('sr'), rgloptions=rgloptions, rglactions=rglactions);

convert -loop 0 fsbrain_subject1_thickness_white.gif fsbrain_subject1_thickness_white_loop_inf.gif

Drawing colorbars

Adding 2D colorbars to a 3D plot has some technical limitations, and fsbrain currently gives you two options to get a colorbar for your plot.

Plotting a colorbar in the background of a 3D plot

The first option is to use the draw_colorbar parameter of the visualization functions that support it (all for which it makes sense e.g., vis.subject.morph.native(), vis.subject.morph.standard(), vis.data.on.subject() and many more). This tries to draw a colorbar into the background of the rgl 3D plot. This is experimental and only works if there is enough space in the plot area. This means that it depends on the views parameter and the size of the rgl device whether or not there is enough space to draw the colorbar. If there is not enough space, the colorbar may not show up at all or it may be empty. You can use the rgloptions parameter to all these functions to increase the size of the rgl device. Example (continued from above):

subjects_dir = fsbrain::get_optional_data_filepath("subjects_dir");
subject_id = 'subject1';
rgloptions=list("windowRect"=c(50, 50, 1000, 1000));     # larger plot
surface = 'white';
measure = 'thickness';
vis.subject.morph.native(subjects_dir, subject_id, measure, 'both', views=c('si'), rgloptions=rgloptions, draw_colorbar = TRUE);

Figure 5: Colorbar integrated into the rgl plot.

Another limitation of this is that the position of the colorbar may not be what you want when using tiled views. It is ok for showing the data to colleagues in a presentation at some informal meeting for most people, but for a figure in a publication, you may want more control over the colorbar appearance and position. In that case, you should use a separate plot for the colorbar.

Plotting a colorbar to a separate 2D plot

This option plots a colorbar for the data of your meshes into a separate 2D plot. You can export the colorbar plot to a graphics file and create the final version of your figure in standard image manipulation software like Inkscape. To do this, you can use the invisible return value of all visualization functions in combination with the coloredmesh.plot.colorbar.separate() function:

subjects_dir = fsbrain::get_optional_data_filepath("subjects_dir");
coloredmeshes = vis.subject.morph.native(subjects_dir, 'subject1', 'thickness', 'lh', views=c('t4'), makecmap_options=list('colFn'=squash::jet));
coloredmesh.plot.colorbar.separate(coloredmeshes, horizontal=TRUE);

Figure 6: Separate horizontal colorbar.

See the help for the coloredmesh.plot.colorbar.separate() function for more customization options.