Downsampling (decimating) a brain surface

Downsampled average cortical surfaces at different iterations (n), with the respective number of vertices (V), edges (E) and faces (F).

In the previous post, a method to display brain surfaces interactively in PDF documents was presented. While the method is already much more efficient than it was when it first appeared some years ago, the display of highly resolved meshes can be computationally intensive, and may make even the most enthusiastic readers give up even opening the file.

If the data being shown has low spatial frequency, an alternative way to display, which preserves generally the amount of information, is to decimate the mesh, downsampling it to a lower resolution. Although in theory this can be done in the native (subject-level) geometry through retessellation (i.e., interpolation of coordinates), the interest in downsampling usually is at the group level, in which case the subjects have all been interpolated to a common grid, which in general is a geodesic sphere produced by subdividing recursively an icosahedron (see this earlier post). If, at each iteration, the vertices and faces are added in a certain order (such as in FreeSurfer‘s fsaverage or in the one generated with the platonic command), downsampling is relatively straightforward, whatever is the type of data.

Vertexwise data

For vertexwise data, downsampling can be based on the fact that vertices are added (appended) in a certain order as the icosahedron is constructed:

  • Vertices 1-12 correspond to n = 0, i.e., no subdivision, or ico0.
  • Vertices 13-42 correspond to the vertices that, once added to the ico0, make it ico1 (first iteration of subdivision, n = 1).
  • Vertices 43-162 correspond to the vertices that, once added to ico1, make it ico2 (second iteration, n = 2).
  • Vertices 163-642, likewise, make ico3.
  • Vertices 643-2562 make ico4.
  • Vertices 2563-10242 make ico5.
  • Vertices 10243-40962 make ico6, etc.

Thus, if the data is vertexwise (also known as curvature, such as cortical thickness or curvature indices proper), the above information is sufficient to downsample the data: to reduce down to an ico3, for instance, all what one needs to do is to pick the vertices 1 through 642, ignoring 643 onwards.

Facewise data

Data stored at each face (triangle) generally correspond to areal quantities, that require mass conservation. For both fsaverage and platonic icosahedrons, the faces are added in a particular order such that, at each iteration of the subdivision, a given face index is replaced in situ for four other faces: one can simply collapse (via sum or average) the data of every four faces into a new one.

Surface geometry

If the objective is to decimate the surface geometry, i.e., the mesh itself, as opposed to quantities assigned to vertices or faces, one can use similar steps:

  1. Select the vertices from the first up to the last vertex of the icosahedron in the level needed.
  2. Iteratively downsample the face indices by selecting first those that are formed by three vertices that were appended for the current iteration, then for those that have two vertices appended in the current iteration, then connecting the remaining three vertices to form a new, larger face.

Applications

Using downsampled data is useful not only to display meshes in PDF documents, but also, some analyses may not require a high resolution as the default mesh (ico7), particularly for processes that vary smoothly across the cortex, such as cortical thickness. Using a lower resolution mesh can be just as informative, while operating at a fraction of the computational cost.

A script

A script that does the tasks above using Matlab/Octave is here: icodown.m. It is also available as part of the areal package described here, which also satisfies all its dependencies. Input and output formats are described here.

Fast surface smoothing on a common grid

Smoothing scalar data on the surface of a high resolution sphere can be a slow process. If the filter is not truncated, the distances between all the vertices or barycentres of faces need to be calculated, in a very time consuming process. If the filter is truncated, the process can be faster, but with resolutions typically used in brain imaging, it can still be very slow, a problem that is amplified if data from many subjects are analysed.

However, if the data for each subject have already been interpolated to a common grid, such as an icosahedron recursively subdivided multiple times (details here), then the distances do not need to be calculated repeatedly for each of them. Doing so just once suffices. Furthermore, the implementation of the filter itself can be made in such a way that the smoothing process can be performed as a single matrix multiplication.

Consider the smoothing defined in Lombardi, (2002), which we used in Winkler et al., (2012):

\tilde{Q}_n = \dfrac{\sum_j Q_j K\left(g\left(\mathbf{x}_n,\mathbf{x}_j\right)\right)}{\sum_j K\left(g\left(\mathbf{x}_n,\mathbf{x}_j\right)\right)}

where \tilde{Q}_n is the smoothed quantity at the vertex or face n, Q_j is the quantity assigned to the j-th vertex or face of a sphere containining J vertices or faces, g\left(\mathbf{x}_n,\mathbf{x}_j\right) is the geodesic distance between vertices or faces with coordinates \mathbf{x}_n and \mathbf{x}_j, and K(g) is the Gaussian filter, defined as a function of the geodesic distance between points.

The above formula requires that all distances between the current vertex or face n and the other points j are calculated, and that this is repeated for each n, in a time consuming process that needs to be repeated over again for every additional subject. If, however, the distances g are known and organised into a distance matrix \mathbf{G}, the filter can take the form of a matrix of the same size, \mathbf{K}, with the values at each row scaled so as to add up to unity, and the same smoothing can proceed as a simple matrix multiplication:

\mathbf{\tilde{Q}} = \mathbf{K}\mathbf{Q}

If the grid is the same for all subjects, which is the typical case when comparisons across subjects are performed, \mathbf{K} can be calculated just once, saved, and reused for each subject.

It should be noted, however, that although running faster, the method requires much more memory. For a filter of full-width at half-maximum (FWHM) f, truncated at a distance t \cdot f from the filter center, in a sphere of radius r, the number of non-zero elements in \mathbf{K} is approximately:

\text{nnz} \approx \dfrac{J^2}{2} \left(1-\cos\left(t \cdot \dfrac{f}{r}\right)\right)

whereas the total number of elements is J^2. Even using sparse matrices, this may require a large amount of memory space. For practical purposes, a filter with width f = 20 mm can be truncated at twice the width (t = 2), for application in a sphere with 100 mm made by 7 subdivisions of an icosahedron, still comfortably in a computer with 16GB of RAM. Wider filters may require more memory to run.

The script smoothdpx, part of the areal interpolation tools, available here, can be used to do both things, that is, smooth the data for any given subject, and also save the filter so that it can be reused with other subjects. To apply a previously saved filter, the rpncalc can be used. These commands require Octave or MATLAB, and if Octave is available, they can be executed directly from the command line.

Figures

The figures above represent facewise data on the surface of a sphere of 100 mm radius, made by recursive subdivision of a regular icosahedron 4 times, constructed with the platonic command (details here), shown without smoothing, and smoothed with filters with FWHM = 7, 14, 21, 28 and 35 mm.

References

Splitting the cortical surface into independent regions

FreeSurfer offers excellent visualisation capabilities with tksurfer and FreeView. However, there are endless other possibilities using various different computer graphics software. In previous posts, it was shown here in the blog how to generate cortical and subcortical surfaces that could be imported into these applications, as well as how to generate models with vertexwise and facewise colours, and even a description of common file formats. It was also previously shown how to arbitrarily change the colours of regions for use with FreeSurfer own tools. However, a method to allow rendering cortical regions with different colours in software such as Blender was missing. This is what this post is about.

The idea is simple: splitting the cortical surface into one mesh per parcellation allows each to be imported as an independent object, and so, it becomes straightforward to apply a different colour for each one. To split, the first step is to convert the FreeSurfer annotation file to a data-per-vertex file (*.dpv). This can be done with the command annot2dpv.

./annot2dpv lh.aparc.annot lh.aparc.annot.dpv

Before running, be sure that ${FREESURFER_HOME}/matlab is in the Octave/matlab, path. With the data-per-vertex file ready, do the splitting of the surface with splitsrf:

./splitsrf lh.white lh.aparc.annot.dpv lh.white_roi

This will create several files names as lh.white_roi*. Each corresponds to one piece of the cortex, in *.srf format. To convert to a format that can be read directly into computer graphics software, see the instructions here.

The annot2dpv and splitsrf are now included in the package for areal analysis, available here.

With the meshes imported, let your imagination and creativity fly. Once produced, labels can be added to the renderings using software such as Inkscape, to produce images as the one above, of the Desikan-Killiany atlas, which illustrates the paper Cortical Thickness or Gray Matter Volume: The Importance of Selecting the Phenotype for Imaging Genetics Studies.

Another method is also possible, without the need to split the cortex, but instead, painting the voxels. This can be done with the command replacedpx, also available from the package above. In this case each region index is replaced by its corresponding statistical value (or any other value), then maps are produced with the dpx2map, shown in an earlier blog post, here. This other method, however, requires that the label indices are known for each region, which in FreeSurfer depends on the rgb colors assigned to them. Moreover, the resulting maps don’t have as sharp and beautiful borders as when the surface is split into independent pieces.

Displaying vertexwise and facewise brain maps

In a previous post, a method to display FreeSurfer cortical regions in arbitrary colours was presented. Suppose that, instead, you would like to display the results from vertexwise or facewise analyses. For vertexwise, these can be shown using tksurfer or Freeview. The same does not apply, however, to facewise data, which, at the time of this writing, is not available in any neuroimaging software. In this article a tool to generate files with facewise or vertexwise data is provided, along with some simple examples.

The dpx2map tool

The tool to generate the maps is dpx2map (right-click to download, then make it executable). Call it without arguments to get usage information. This tool uses Octave as the backend, and it assumes that it is installed in its usual location (/usr/bin/octave). It is also possible to run it from inside Octave or Matlab using a slight variant, dpx2map.m (in which case, type help dpx2map for usage).

In either case, the commands srfread, dpxread and mtlwrite must be available. These are part of the areal package discussed here. And yes, dpx2map is now included in the latest release of the package too.

To use dpx2map, you need to specify a surface object that will provide the geometry on which the data colours will be overlaid, and the data itself. The surface should be in FreeSurfer format (*.asc or *.srf), and the data should be in FreeSurfer “curvature” format (*.asc, *.dpv) for vertexwise, or in facewise format (*.dpf). A description of these formats is available here.

It is possible to specify the data range that will be used when computing the scaling to make the colours, as well which range will be actually shown. It is also possible to split the scale so that a central part of it is omitted or shown in a colour outside the colourscale. This is useful to show thresholded positive and negative maps.

The output is saved either in Stanford Polygon (*.ply) for vertexwise, or in Wavefront Object (*.obj + *.mtl) for facewise data, and can be imported directly in many computer graphics software. All input and output files must be/are in their respective ascii versions, not binary. The command also outputs a image with the colourbar, in Portable Network Graphics format (*.png).

An example object

With a simple geometric shape as this it is much simpler to demonstrate how to generate the maps, than using a complicated object as the brain. The strategy for colouring remains the same. For the next examples, an ellipsoid was created using the platonic command. The command line used was:

platonic ellipsoid.obj ico sph 7 '[.25 0 0 0; 0 3 0 0; 0 0 .25 0; 0 0 0 1]'

This ellipsoid has maximum y-coordinate equal to 3, and a maximum x- and z-coordinates equal to 0.25. This file was converted from Wavefront *.obj to FreeSurfer ascii, and scalar fields simply describing the coordinates (x,y,z), were created with:

obj2srf ellipsoid.obj > ellipsoid.srf
srf2area ellipsoid.srf ellipsoid.dpv dpv
gawk '{print $1,$2,$3,$4,$2}' ellipsoid.dpv > ellipsoid-x.dpv
gawk '{print $1,$2,$3,$4,$3}' ellipsoid.dpv > ellipsoid-y.dpv
gawk '{print $1,$2,$3,$4,$4}' ellipsoid.dpv > ellipsoid-z.dpv

It is the ellipsoid-y.dpv that is used for the next examples.

Vertexwise examples

The examples below use the same surface (*.srf) and the same curvature, data-per-vertex file (*.dpv). The only differences are the way as the map is generated and presented, using different colour maps and different scaling. The jet colour map is the same available in Octave and Matlab. The coolhot5 is a custom colour map that will be made available, along with a few others, in another article to be posted soon.

Example A

In this example, defaults are used. The input files are specified, along with a prefix (exA) to be used to name the output files.

dpx2map ellipsoid-y.dpv ellipsoid.srf exA

Example B

In this example, the data between values -1.5 and 1.5 is coloured, and the remaining receive the colours of the extreme points (dark blue and dark red).

dpx2map ellipsoid-y.dpv ellipsoid.srf exB jet '[-1.5 1.5]'

Example C

In this example, the data between -2 and 2 is used to define the colours, with the values below/above receiving the extreme colours. However, the range between -1 and 1 is not shown or used for the colour scaling. This is because the dual option is set as true as well as the coption.

dpx2map ellipsoid-y.dpv ellipsoid.srf exC coolhot5 '[-2 2]' '[-1 1]' true '[.75 .75 .75]' true

Example D

This example is similar as above, except that the values between -1 and 1, despite not being shown, are used for the scaling of the colours. This is due to the coption being set as true.

dpx2map ellipsoid-y.dpv ellipsoid.srf exD coolhot5 '[-2 2]' '[-1 1]' true '[.75 .75 .75]' false

Example E

Here the data between -2 and 2 is used for scaling, but only the points between -1 and 1 are shown. This is because the option dual was set as false. The values below -1 or above 1 receive the same colours as these numbers, because the coption was configured as true. Note that because all points will receive some colour, it is not necessary to define the colourgap.

dpx2map ellipsoid-y.dpv ellipsoid.srf exE coolhot5 '[-2 2]' '[-1 1]' false '[]' true

Example F

This is similar as the previous example, except that the values between -1 and 1 receive a colour off of the colour map. This is because both dual and coption were set as false.

dpx2map ellipsoid-y.dpv ellipsoid.srf exF coolhot5 '[-2 2]' '[-1 1]' false '[.75 .75 .75]' false

Facewise data

The process to display facewise data is virtually identical. The only two differences are that (1) instead of supplying a *.dpv file, a *.dpf file is given to the script as input, and (2) the output isn’t a *.ply file, but instead a pair of files *.obj + *.mtl. Note that very few software can handle thousands of colours per object in the case of facewise data. Blender is recommended over most commercial products specially for this reason (and of course, it is free, as in freedom).

Download

The dpx2map is available here, and it is also included in the areal package, described here, where all its dependencies are satisfied. You must have Octave (free) or Matlab available to use this tool.

How to cite

If you use dpx2map for your scientific research, please, remember to mention the brainder.org website in your paper.

Update: Display in PDF documents

3D models as these, with vertexwise colours, can be shown in interactive PDF documents. Details here.

Merging multiple surfaces

Say you have a number of meshes in FreeSurfer ascii format (with extension *.asc or *.srf), one brain structure per file. However, for later processing or to import in some computer graphics software, you would like to have these multiple meshes all in a single file. This post provides a small script to accomplish this: mergesrf.

To use it, right click and save the file above, make it executable and, ideally, put it in a place where it can be found (or add its location to the environmental variable ${PATH}. Then run something as:

mergesrf file1.srf file2.srf fileN.srf mergedfile.srf

In this example, the output file is saved as mergedfile.srf. Another example is to convert all subcortical structures into just one large object, after aseg2srf as described here. To convert all, just change the current directory to ${SUBJECTS_DIR}//ascii, then run:

mergesrf * aseg_all.srf

A list with the input files and the output at the end is shown below:

The script uses Octave, which can be downloaded freely. The same script, with a small modification, can also run from inside matlab. This other version can be downloaded here: mergesrf.m

Requirements

In addition to Octave (or matlab), the script also requires functions to read and write surface files, which are available from the areal package (described here and downloadable here).

Brain meshes available

A set of 3D brain meshes produced with FreeSurfer and individually partitioned into separate files following the atlases of Desikan-Killiany, Destrieux et al., and Desikan-Killiany-Tourville (DKT), is now available for download here. Surfaces for subcortical structures are also available.

These meshes are meant to be used to help with scientific visualisation and/or artistic rendering in computer graphics software, most prominently Blender, but also in any other application that can import files in Wavefront (*.obj) or Stanford Polygon (*.ply) formats. The released files are under a Creative Commons license. See the download page for details.

Facewise brain cortical surface area: scripts now available

https://brainder.org/download/areal

In the paper Measuring and comparing brain cortical surface area and other areal quantities (Neuroimage, 2012;61(4):1428-43, doi:10.1016/j.neuroimage.2012.03.026, pmid:22446492), we described the steps to perform pycnophylactic interpolation, i.e., mass-conservative interpolation, of brain surface area. The method is also valid for any other quantity that is also areal by nature.

The scripts to perform the steps we propose in the paper are now freely available for Octave and/or matlab. They can also be executed directly from the shell using Octave as the interpreter. In this case, multiple instances can be easily submitted to run in a cluster. Click here to dowload and for step-by-step instructions.

Importing FreeSurfer subcortical structures into Blender

In a previous article, a method to import cortical meshes from FreeSurfer into Blender and other CG applications was presented. Here it is complemented with a method to import the subcortical structures. These structures are segmented as part of FreeSurfer subcortical stream, and are stored in a volume-based representation. In order to import into Blender, first it is necessary to generate surfaces, then convert to Wavefront OBJ, and finally, do the import. Below is an example of how the output looks like.

An example of a subcortical meshes generated from FreeSurfer subcortical segmentations imported into Blender.

Generating the surfaces for subcortical structures

The surfaces can be generated using FreeSurfer’s own tools. The script aseg2srf automates the task. Suppose you would like to create meshes for all subcortical structures of the subjects named “capitu”, “bentinho” and “ezequiel”. Simply run:

aseg2srf -s "capitu bentinho ezequiel"

The outputs are saved in a directory called ascii, inside each subject’s directory. The script requires that FreeSurfer is installed and working properly, and that the environmental variable ${SUBJECTS_DIR} has been set to point to the location where the subjects are. It is also necessary to have writing permissions to this directory, so that the surfaces can be saved.

If there are many subjects, provide a list from a text file, with one subject per line:

aseg2srf -s list.txt

To make surfaces for just some structures, provide their label numbers with the -l option. For the hippocampus, for instance, the labels are 17 and 53 for left and right respectively:

aseg2srf -s "capitu bentinho ezequiel" -l "17 53"

The script has also a -d option, for debugging, which allows to keep temporary files used during the mesh generation.

List of labels

Below is a list of the label codes for the major structures. A more comprehensive list is in the aseg.stats file, inside the stats directory for each subject, or in the FreeSurferColorLUT.txt file (in the FreeSurfer directory).

Left lateral ventricle 4
Left inferior part of the lateral ventricle 5
Left cerebellum white matter 7
Left cerebellum cortex 8
Left thalamus 10
Left caudate 11
Left putamen 12
Left pallidum 13
Left hippocampus 17
Left amygdala 18
CSF 24
Left accumbens 26
Left ventral diencephalon 28
Left choroid plexus 31
Right lateral ventricle 43
Right inferior part of the lateral ventricle 44
Right cerebellum white matter 46
Right cerebellum cortex 47
Right thalamus 49
Right caudate 50
Right putamen 51
Right pallidum 52
Right hippocampus 53
Right amygdala 54
Right accumbens 58
Right ventral diencephalon 60
Right choroid plexus 63
Corpus callosum, posterior part 251
Corpus callosum, middle posterior part 252
Corpus callosum, central part 253
Corpus callosum, middle anterior part 254
Corpus callosum, anterior part 255
3rd ventricle 14
4th ventricle 15
5th ventricle (CSF sometimes found in the septum pellucidum) 72
Brain stem 16

Conversion to Wavefront OBJ

Once the meshes have been generated, the process of conversion into a format that can be read by CG software is the same as described previously here.

Download

The script is available here: aseg2srf.

Requirements

FreeSurfer must be installed and working properly. The environmental variable ${SUBJECTS_DIR} must point to the directory where the data for each subject is located. It is also assumed that recon-all has been run successfully.

Importing FreeSurfer cortical meshes into Blender

Computer graphics (CG) applications, such as Blender, Maya or 3D Studio, doesn’t come with filters that allow to directly import or export FreeSurfer meshes. To import brain surfaces to these applications, converting first to an intermediate file format solves the problem. Below is an example of a brain imported into Blender.

An example of a brain mesh from FreeSurfer imported into Blender after conversion to Wavefront OBJ.

The mris_convert, from FreeSurfer, allows conversion to stereolithography format (STL). However, not all CG applications read natively from STL files. Moreover, other formats, such as Wavefront OBJ and Stanford PLY have simpler structure — yet allowing features that will be discussed here in forthcoming articles.

Conversion from FreeSurfer to Wavefront OBJ

A simple procedure to convert a FreeSurfer surface to Wavefront OBJ (ASCII) is:

1) Use mris_convert to convert from FreeSurfer binary to ASCII:

mris_convert lh.pial lh.pial.asc

2) A suggestion is to rename the .asc file to .srf, to differentiate it from other files that also have .asc extension and have a different internal structure.

mv lh.pial.asc lh.pial.srf

3) Use the script srf2obj to convert from .srf to .asc:

srf2obj lh.pial.srf > lh.pial.obj

The symbol “>” is necessary so that the output from srf2obj is redirected to a new file.

Conversion from Wavefront OBJ to FreeSurfer

To convert back to FreeSurfer, a similar strategy is used.

1) Use the script obj2srf to generate new .srf files (here directly named as .asc, as expected by mris_convert)

obj2srf lh.pial.obj > lh.pial.asc

2) Use mris_convert:

mris_convert lh.pial.asc lh.pial

Download

Both scripts are available here: srf2obj and obj2srf. Right-click and save, then make them executable.

Requirements

The scripts use gawk as the interpreter, and expects to find it on its usual location (/bin/gawk). If in your system gawk is installed elsewhere, modify the first line accordingly.

If you don’t have gawk, you can get it from the repository of your Linux distribution. For Mac, gawk is available in MacPorts and Fink. Alternatively, you may replace gawk with awk itself, nawk, or any other awk implementation that works.

A script to create platonic polyhedra, with recursive subdivision

Suppose you need a tri-dimensional file that contains a platonic polyhedron, say, an octahedron or an icosahedron, stored as a mesh. This file should, however, meet certain criteria: be in an useful format that can be imported into modelling applications, such as Blender, and be easily manipulable using scripts in languages suited for numerical computations, such as Octave or R and perhaps even be manipulated directly from the shell, using, e.g., gawk.

To meet these demands, below is a tool that creates any of the five platonic polyhedra in Wavefront Object format, with the *.obj extension (in ASCII). The user can specify the desired edge length, the area of the face, the total area, the total volume, the radius of the circumscribed sphere, or the radius of the inscribed sphere.

The tool also allows recursive subdivision of the faces of the polyhedra that are made of triangular faces (tetrahedron, octahedron and icosahedron). The subdivision implemented allows geodesic spheres of Class I (Kenner, 1976) to be produced easily. It is also possible to apply affine transformations to the coordinates of the vertices, so that the resulting mesh can be scaled, translated, rotated and sheared in any direction. Finally, it is possible to randomly perturb the vertex positions, so that the mesh becomes irregular, yet preserving the original topology.

In the case of recursive subdivision, it is important to note that not all edges have the same length, and not all faces have the same area. See the figure below for an illustration:

(a) A geodesic sphere can be produced from recursive subdivision of a regular icosahedron. At each iteration, the number of faces is quadruplied. (b) After the first iteration the faces no longer have regular sizes, with the largest face being approximately 1.3 times larger than the smallest as n increases.

The formulas to calculate the number of vertices, edges and faces are:

  • Number of vertices: V=4^n(V_0-2)+2;
  • Number of edges: E=4^nE_0;
  • Number of faces: F=4^nF_0;

where F_0, V_0 and E_0 are, respectively, the number of faces, vertices and edges of the polyhedron with triangular faces used for the initial subdivision:

  • Tetrahedron: V_0=4, E_0=6 and F_0=4;
  • Octahedron: V_0=6, E_0=12 and F_0=8;
  • Icosahedron: V_0=12, E_0=30 and F_0=20.

The figure above and the formulas have been discussed by Winkler et al., 2012.

Main file

The links to the script are below. Either of these will work; choose according to your needs:

  • platonic – This is the Octave script. Simply download it, change the first line to point to your correct Octave location, make it executable, and run it directly from the shell. Call it without arguments to obtain usage information.
  • platonic.m – This is the same as above, but can be executed from within Octave or MATLAB. Type ‘help platonic’ to obtain usage information.

Requirements

  • subdivtri.m – Add this to your Octave or MATLAB path.

Examples

To create an icosahedron with radius 100, recursively subdivided 7 times, if from the shell, use:

platonic ico7.obj ico sph 7 100

or, if from inside Octave/MATLAB, use:

platonic('ico7.obj','ico','sph',7,100);

The outputs from platonic are always in Wavefront (.obj) format. To use this icosahedron, e.g., with areal interpolation, it needs to be converted to FreeSurfer ascii format:

obj2srf ico7.obj > ico7.srf

The obj2srf command is available in the package for areal analysis, here.

References