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

Quickly inspect FreeSurfer cortical surfaces

To quickly inspect the surface reconstruction and parcellations produced by the cortical stream of FreeSurfer, rather than laboriously open each hemisphere at a time in tksurfer, an easier solution is to run a script that captures images, organise and present them all in a single html page that can be seen in any browser.

Before you start

Before you start, make sure you have:

  • FreeSurfer, any version released in the last 2 years. The procedure has been tested up to version 5.1.0.
  • ImageMagick, any recent version. For most Linux, ImageMagick can be installed directly from the distro’s repository. For Mac, it can be obtained with MacPorts or Fink.
  • The scripts takeshots, makehtml, shots_tksurfer.tcl, and shots_tkmedit.tcl. Put all these scripts into the same directory and make takeshots and makehtml executable.

Capturing the images

  1. Make sure recon-all finished for all your subjects. If there are no results, obviously there won’t be anything to inspect.
  2. Make sure the variable SUBJECTS_DIR is correctly set, pointing to the directory that contains the subjects, and that you have writing permissions to this directory.
  3. Prepare a text file containing a list of subjects, one per line, like this:
    subj001
    subj002
    ...
    subjNN
  4. Run the takeshots script. A typical usage is like this:
    ./takeshots -l listsubj.txt -m pial -m inflated -p aparc
    The options are:
    -l: List of subjects.
    -s: Additional subjects not in the list.
    -m: Surfaces to inspect (pial, inflated, white, sphere, etc).
    -p: Parcellations to inspect (aparc, aparc.a2005s, aparc.a2009s, etc).
    This script will create a subdirectory called ‘shots’ inside the directory of each subject, and will store a number of .tif files that are the different views of the cortical surfaces.
  5. Run the makehtml script. A typical usage is like this:
    ./makehtml -l listsubj.txt -m pial -m inflated -p aparc -d /path/to/html
    The options are identical to the previous command, with the addition of the option -d that lets you specify a directory where the html will be created. If this directory doesn’t exist, it will be created.

What you get

Open the file index.html in the directory specified in the previous step. There you’ll see all the surfaces of all the subjects, all in the same page, allowing you to quickly compare and have an idea of how good or bad they are. To see the name of the subject, hover the mouse on the corresponding image. Click on the image to see it larger. An example, using only pial and inflated surfaces, and the aparc parcellation, is here.

FreeSurfer brains in arbitrary colours

Suppose you run a statistical test for each of the regions of the parcellated brain produced by FreeSurfer. It can be either surface area comparison between groups, or maybe differences in thickness. For each region, you may obtain a statistical score and likely a p-value as well. How to display these results in a color-coded brain model, still inside FreeSurfer?

The process is consists of three steps:

  1. Create of a table with your results, in a format that FreeSurfer can read.
  2. Embed the table in an annotation file.
  3. Display in tksurfer.

1. Create a table

The table must contain not the actual scalars, but instead rgb triplets. This can be done using a simple Octave or matlab script. Since often these results are tabulated as spreadsheets, an alternative, straightforward way is using software like LibreOffice, which gives instantaneous graphical feedback of how the resulting table looks like.

An example of such a spreadsheet is this: colormap.ods. The data go into the sheets named lh.data and rh.data. For structures that are not supposed to display any color, put an off-range value below the minimum value of the other regions (e.g. -1 if the scalars are all positive). Do not add or remove structures.

The sheet named colormap contains the key between the scalar values and the actual rgb colors. In this sheet, the first column (column A) is simply a rank number, used in the calculations; the second column (col. B) contains the range of values between 0 and max, regularly spaced. The remaining three columns (C, D and E) contains the rgb triplets to be assigned to each value in column B. For each scalar to be shown, the closest number has its corresponding rgb used.

The sheets named lh.aparc.annot.ctab and rh.aparc.annot.ctab contain the resulting tables that will be used for actual display purposes. Each contain an index number, starting from 0, the structure name as used in FreeSurfer, the rgb colors, and a fourth value that is called “flag” and is usually zero.

In most cases, only the numerical values in lh.data and rh.data have to be changed. In some occasions, it may be necessary to change the range to be shown (i.e., column B on the colormap sheet), or the colormap itself (columns C, D and E). Different colormaps can be generated, for instance, using Octave. In the file used here, the colormap is the “spectrum”, also known as “jet”, shown below:

Finally, FreeSurfer cannot read binary spreadsheets like this. The file has to be converted to simple values separated by space. In LibreOffice or OpenOffice, there is a very simple way of doing this. Go to File -> Save As… In the File type, choose Text CSV, and make sure to mark the option Edit Filter Settings:

On the next window, as Field Delimiter, choose {space}. As Text delimiter, leave blank, as shown below:

The resulting file, saved to the disk, should look like this (click here for the full file):

0 bankssts 92 255 163 0
1 caudalanteriorcingulate 143 255 112 0
2 caudalmiddlefrontal 0 92 255 0
[...]
34 transversetemporal 255 133 0 0
35 unknown 160 160 160 0

2. Embed the table into annotation file

Once a table has been produced, it can be embedded into the annotation file that defines the the vertices that belong to each region. The annotation files produced by FreeSurfer already contain a default color table, which can be replaced. This can be accomplished by the Octave function replace_ctab.m. This function needs Octave 3.4.0+ or matlab. Type help replace_ctab at the prompt to get usage information. Here is an example:

replace_ctab('${SUBJECTS_DIR}/bert/label/lh.aparc.annot',...
'${SUBJECTS_DIR}/bert/label/lh.aparc.annot.ctab.csv',...
'${SUBJECTS_DIR}/bert/label/lh.aparc.annot.new');

In this example, the original annotation is stored in the file ${SUBJECTS_DIR}/bert/label/lh.aparc.annot, the new color table, as saved in LibreOffice, is in the file ${SUBJECTS_DIR}/bert/label/lh.aparc.annot.ctab. The new annotation file, with the new colortable, will be in the file ${SUBJECTS_DIR}/bert/label/lh.aparc.annot.new. Notice that for the function replace_ctab to work, the directory ${FREESURFER_HOME}/matlab must be in the Octave or matlab path.

3. Display with tksurfer

Once the new annotation file has been produced, display the result is straightforward:

tksurfer bert lh pial -annotation lh.aparc.annot.new

The result should be like the brain shown at the beginning of this article.

Update (01.Aug.2013): In more recent FS versions, the “lh” in front of the annotation file can be omitted, as it’s already specified in the same command line (thanks to Krishna Pancholi, Olin Neuropsychiatric Research Center, for the tip.). This affects only the last command line, which becomes then:

tksurfer bert lh pial -annotation aparc.annot.new