Here a software implementation of the steps to perform areal (pycnophylactic) interpolation of brain surface area is available, using Octave/matlab. The method can be used for area itself, or for of any other quantity that needs mass-conservative interpolation, including cortical volumes and possibly other measurements. These steps have been published in:
- Winkler AM, Sabuncu MR, Yeo BT, Fischl B, Greve DN, Kochunov P, Nichols TE, Blangero J, Glahn DC. Measuring and comparing brain cortical surface area and other areal quantities. Neuroimage. 2012 Jul 16;61(4):1428-43.
To get the most recent version of the scripts, visit the “areal” repository in GitHub.
There are two directories in the package:
share: This contains some functions (m-files) that run from inside Octave or matlab. Simply add this directory to your
~/.octavercfile or to matlab‘s path. Type
helpat the prompt obtain usage information.
bin: This contains copies of most of the same functions, but as executable scripts that run directly from the shell, using Octave as the interpreter. You may consider adding this directory to your
sharedirectory still needs to be in the
~/.octavercfile. To obtain usage information, call each command without arguments.
If you run using Octave, make sure to use version 3.4.0 or more recent, and that the command
convhull is available (i.e. Octave was compiled with the Qhull library). Moreover, the scripts assume that Octave executable is on its usual location, i.e.,
/usr/bin/octave. If you have it installed elsewhere, edit the first line of the scripts accordingly (or make a symbolic link if you have administrative privileges in your system).
- Run FreeSurfer in all the subjects that will enter into the statistical analysis;
- Run Spherical Demons for registration to a common space;
- Compute the area-per-face from native geometry;
- Perform the areal interpolation;
- (Optional) Convert to vertexwise;
- Run the statistical analysis.
Step 1: FreeSurfer
It is not mandatory to use FreeSurfer to construct the surfaces, and other tools are available (e.g. Caret, BrainVISA, civet and others). However, here it is assumed that the data come from FreeSurfer. If you choose to use other software packages to create surfaces, make the appropriate file conversions are done as needed so that the subsequent steps run correctly.
Step 2: Spherical Demons
Follow the instructions in the Spherical Demons website. SD will create a subdirectory called
SD inside each subject’s directory, and will also create files named
?h.sphere.SD.reg in the
surf subdirectory for each subject. These registered spheres will be used during interpolation.
Step 3: Facewise area
The tools work with ascii file formats. The reason is that it easier to identify eventual errors, to perform quality control, and to transfer between different applications. A description of the file formats is here. The script to perform batch conversion between file formats,
subj2ascii, is included in the package above.
Once all relevant surface and curvature files for all subjects have been converted, it’s time to compute the facewise surface area. This is done with the command
srf2area. From inside Octave or matlab, for each subject, use something as:
srf2area('lh.white.srf', 'lh.white.dpf', 'dpf');
lh.white.srf is the surface file in ascii format, and
lh.white.dpf is the file that will be created. The
dpf option means that a data-per-face file should be created. If running directly from the shell, use
srf2area lh.white.srf lh.white.dpf dpf
You may consider using a for-loop, either inside Octave, or directly in the shell. The calculation of the facewise area is generally fast.
Step 4: Areal interpolation
The command to run the interpolation is
srf2srf. It contains 4 different interpolation methods.
areal: This is probably the one you want, and it is the same we described in the paper. It is the recommended method for analysis of brain surface area. The areal interpolation is pycnophylactic, i.e., it conserves the areal quantities locally, regionally and globally. The amount of area each target face receives is determined by the overlap with each of the source faces. See the paper for detailed description.
barycentric: This will perform simple barycentric interpolation. It is not a mass-conservative method and so, cannot be used for brain surface area. It is nonetheless available as it can be used for thickness, and also as an intermediate step if one wants to retessellate a surface to a different geometry.
distributive: This is a vertexwise mass-conservative method that works by splitting the area in each vertex using the proportion given by the barycentric coordinates of the source vertex in relation to the target face, and redistributing these quantities across the three vertices in the target, adding to what was already there. At the beginning of the loop, each target vertex has zero areal quantity, and at the end, each will have received an amount that is taken from the nearby vertices from the sources using the proportion as calculated above. It isn’t a bad method, but the
arealis far superior, hence it is the one we advocate and showed in the paper.
nearest: This is a simple nearest-neighbour interpolation. It was used only for testing purposes and left here if anyone needs. However, if you really need nearest-neighbour interpolation, use instead the one that already comes with FreeSurfer, because that one uses a compiled routine that runs much faster.
To run areal interpolation from inside Octave or matlab use:
srf2srf('areal', 'lh.sphere.SD.reg.srf', 'ico7.srf', 'lh.overlaps.csv', 'lh.white.dpf', 'lh.white.ico7.dpf');
areal is the interpolation method,
lh.sphere.SD.reg.srf is the source surface,
ico7.srf is the target surface,
lh.overlaps.csv is a table that will be created containing the areas of overlap between each source and target face,
lh.white.dpf is the source areal quantity, here the facewise area from native geometry as calculated in the previous step, and
lh.white.ico7.dpf is the output facewise area after interpolation, i.e., in the resolution of the
ico7.srf. If running directly from the shell, use:
srf2srf areal lh.sphere.SD.reg.srf ico7.srf lh.overlaps.csv lh.white.dpf lh.white.ico7.dpf
ico7.srf can be created using the
platonic command, as described here. Alternatively, use either
rh.sphere.reg (these are the same) from the FreeSurfer
fsaverage. The benefit of using these files is that the results can be more easily projected in the
fsaverage for visualisation in FreeSurfer, although during to small rounding errors, the fs spheres are noisier. If you choose to use the
ico7.srf created with
platonic, retessellated versions of
fsaverage are available here.
The resulting table
lh.overlaps.csv can be used to interpolate quickly other areal quantities from the same source and to the same target surfaces, without having to re-run
srf2srf again, as the last can be slow, particularly in old hardware. The command to accomplish this task is
applyolp, also inside the download package.
Step 5: Smooth
To smooth the interpolated file, use the command
smoothdpx. However, before smoothing, it’s necessary to consider that, even for the subdivided icosahedron, the faces have not all the same size. To account for this, use the simple matlab calculator provided in the package. From the shell, the call is:
rpncalc lh.white.ico7.dpf load ico7.dpf load '/' 125663.706144 '*' 327680 '/' lh.white.ico7.tmp.dpf save
or, from inside Octave/matlab:
rpncalc('lh.white.ico7.dpf', 'load', 'ico7.dpf', 'load', '/', '125663.706144','*', '327680', '/', 'lh.white.ico7.tmp.dpf', 'save');
Note that, even from inside matlab, the numeric arguments need to be provided as strings, i.e., between single quotes. The reason is to allow compatibility with the command-line version. The ‘ico7.dpf’ is the area-per-face of the ico7, which can be calculated with the command
srf2area The number 125663.706144 is , and 327680 is the number of faces in the geodesic sphere (see the Appendix a of the paper). Once the different face sizes have been taken into account, smoothing can be performed. From Octave or matlab, prompt, the syntax is:
smoothdpx('-i', 'lh.white.ico7.tmp.dpf', '-s', 'ico7.srf', '-o', 'lh.white.ico7.fwhm10.dpf', '-f', 10);
where 10 represents the fwhm in mm. If calling directly from the shell, the syntax is:
smoothdpx -i lh.white.ico7.tmp.dpf -s ico7.srf -o lh.white.ico7.fwhm10.dpf -f 10
Smoothing uses the moving weights method, and the ico7 surface file is necessary so that the command knows the geometry and information about the neighbourhood between faces. The command
smoothdpx can save the smoothing filter in matrix form, which can be used to smooth quickly other subjects that have resampled to the same common grid; the smoothing is then applied with matrix multiplication using
For certain statistical analysis, after smoothing you may consider to undo the mathematical operation done just before smoothing. The
rpncalc can be used again for that.
Step 6: (Optional) Convert to vertexwise
If only software that can display vertexwise data is available, or if the study will include data that is vertexwise (e.g. thickness, fMRI projected to the surface), then conversion from facewise to vertexwise is necessary. The conversion needs to be done after areal interpolation, but before statistical analysis.
The command to convert, from the shell, is:
dpf2dpv ico7.srf lh.white.ico7.dpf lh.white.ico7.dpv
or, from inside Octave/matlab:
dpf2dpv('ico7.srf', 'lh.white.ico7.dpf', 'lh.white.ico7.dpv');
It is immaterial whether it is performed before or after smoothing. Since the conversion halves the resolution, run it before smoothing makes the smoothing faster. However, it is also then necessary to convert the face sizes of the ico7 (ico7.dpf) to vertexwise. The
rpncalc command works in the same way for both vertexwise or facewise data.
Step 7: Statistical analysis
Use your favourite software for the statistical analysis. If you followed the above steps, at this point all files are in ascii format, hence load data into statistical packages should not be an issue. Each software has its own method to load data.
The Box–Cox transformation, that is discussed in the paper, is available as a command on its own, e.g., in matlab, r, spss, as well as in other packages.
- 02.Jun.2012: First public release.
- 23.Aug.2012: Included the tool
rpncalcto take different face sizes into account, as described in the Appendix a of the paper.
- 25.Aug.2012: Fixed reading of curvature files (
srf2srf); fixed use of non-char values (
platonic). Thanks to Gregory Kirk (University of Wisconsin–Madison).
- 26.Aug.2012: Added a missing dependency for
- 28.Aug.2012: Added step for conversion from facewise to vertexwise.
- 21.Sep.2012: Fixed crashes with
convhullin some matlab versions. Thanks to Inge Amlien (Universitetet i Oslo).
- 24.Sep.2012: Added a missing dependency for the smoothing step. Fixed reading of curvature files (
dpx2dpv). Thanks to Gregory Kirk (University of Wisconsin–Madison).
- 02.Oct.2012: Fixed an issue in
platonicthat affected interpolation to spheres of small radii. Thanks to Christopher Owen (Washington University in St Louis).
- 09.Jun.2013: Fixed call to
srf2area. Thanks to Knut Bjuland (Norwegian University of Science and Technology).
- 28.Jul.2013: Included a tool to display the results,
dpx2map, and its dependencies. This tool is described here.
- 08.Jan.2014: Included tools to split the surface and to replace color indices, as described here.
- 03.Feb.2014: Fixed issue with argumernt parsing in
- 25.Jul.2014: Added the retessellated
fsaveragesurfaces for use with
ico7. Thanks to André Zugman (King’s College London), the first to bump into this problem.
- 28.Nov.2014: Updated
rpncalcwith more functions, for much faster smoothing, as described here.
- 05.Jul.2015: Repository now on GitHub.
- 31.May.2016: Added tools
ply2idtfto downsample the icosahedron and insertion into PDF documents (details here and here).