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.

The paper is available on PubMedCentral (Open Access) or on ScienceDirect.

Download

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 ~/.octaverc file or to matlab‘s path. Type help <function name> at 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 PATH variable. The share directory still needs to be in the ~/.octaverc file. 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).

Instructions

1. Run FreeSurfer in all the subjects that will enter into the statistical analysis;
2. Run Spherical Demons for registration to a common space;
3. Compute the area-per-face from native geometry;
4. Perform the areal interpolation;
5. Smooth;
6. (Optional) Convert to vertexwise;
7. Run the statistical analysis.

Step 1: FreeSurfer

Run recon-all on all the subjects that will be analysed. Inspect the resulting surfaces to make sure that the result is acceptable. You can use the QAtools, or go directly to the surfaces.

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');

where 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 areal is 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');

where 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

The ico7.srf can be created using the platonic command, as described here. Alternatively, use either lh.sphere.reg or 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 rpncalc.

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.

Version history

• 02.Jun.2012: First public release.
• 23.Aug.2012: Included the tool rpncalc to 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 platonic.
• 28.Aug.2012: Added step for conversion from facewise to vertexwise.
• 21.Sep.2012: Fixed crashes with qhull/convhull in 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 platonic that affected interpolation to spheres of small radii. Thanks to Christopher Owen (Washington University in St Louis).
• 09.Jun.2013: Fixed call to dpxread in 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 applyolp.
• 25.Jul.2014: Added the retessellated fsaverage surfaces for use with platonic‘s ico7. Thanks to André Zugman (King’s College London), the first to bump into this problem.
• 28.Nov.2014: Updated smoothdpx and rpncalc with more functions, for much faster smoothing, as described here.
• 05.Jul.2015: Repository now on GitHub.
• 31.May.2016: Added tools icodown and ply2idtf to downsample the icosahedron and insertion into PDF documents (details here and here).