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. Typehelp
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 yourPATH
variable. Theshare
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
- 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;
- Smooth;
- (Optional) Convert to vertexwise;
- 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 theareal
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
insrf2area
. 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 withplatonic
‘sico7
. Thanks to André Zugman (King’s College London), the first to bump into this problem. - 28.Nov.2014: Updated
smoothdpx
andrpncalc
with more functions, for much faster smoothing, as described here. - 05.Jul.2015: Repository now on GitHub.
- 31.May.2016: Added tools
icodown
andply2idtf
to downsample the icosahedron and insertion into PDF documents (details here and here).