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

Braindering with ASCII files

Binary file formats are faster to read and occupy less space for storage. However, when developing tools or when analysing data using non-standard methods, ascii files are easier to treat, and allow certain errors to be spotted immediately. Working with ascii files offer a number of advantages in these circumstances: the content of the files are human readable, meaning that there is no need for any special program other than perhaps a simple text editor, and can be manipulated easily with tools as sed or awk. It should be noted though, that while binary files require attention with respect to endianess and need some metadata to be correctly read, ascii need proper attention to end-of-line characters, that vary according to the operating system.

In articles to be posted soon I hope to share some tools to analyse and visualise brain imaging data. To use these tools, the files must be in ascii format. Each of these formats is discussed below and an example is given for a simple octahedron.

FreeSurfer surface (*.asc | suggested: *.srf)

FreeSurfer natively saves its surface files in binary form. The command mris_convert allows conversion to ascii format. The file extension is *.asc. This extension, however, is ambiguous with the curvature format (see below), reason for which a recommendation is to rename from *.asc to *.srf immediately after conversion.

The *.asc/*.srf format has a structure as described below. This format does not allow storing normals, colors, materials, multiple objects or textures. Just simple geometry.

FreeSurfer curvature (*.asc | suggested: *.dpv)

A curvature file contains a scalar value assigned to each vertex. This file format does not store the geometry of the surface and, in binary form, a curvature file simply contains one scalar per vertex. In ascii format, the vertex coordinates are also saved, but not how these vertices relate to define faces and the overall geometry.

Like with the FreeSurfer surface files, curvature files in ascii format also have the extension *.asc. To avoid confusion, it is suggested to rename from *.asc to *.dpv (data-per-vertex) immediately after conversion.

Facewise data (*.dpf)

This is not a native format from FreeSurfer, but a modification over the ascii curvature file described above. The file stores a scalar value per face, i.e. data-per-face (dpf) and, instead of saving vertex coordinates as with the *.dpv files, here the vertex indices that define the faces are saved. It is done like so because if both *.dpv and *.dpf files are present for a given geometric object, the object can be reconstructed with very little manipulation using simply awk. Currently there is no binary equivalent format.

VTK polygonal data (*.vtk)

A number of file formats have been developed over the years for use with the Visualization Toolkit (vtk). Some of these formats are now considered as “legacy” and are being phased out in favour of xml formats. Despite this, they are still used by a number of applications, such as fsl/first. The octahedron above in vtk polygonal data (polydata) format would have a structure as depicted below.

Wavefront Object (*.obj)

Wavefront Object files have a very simple structure, yet are powerful to allow definition of textures, normals and different materials. The format is supported by a variety of different graphic applications. The octahedron above would have a file structure as shown below.

The example does not show colors. However, it is possible to assign a material to each face, and materials contain their own color. These materials are typically stored in a separate file, a Materials Library, with extension *.mtl. An example for the same octahedron, both the object file and its associated materials library, is here: octa-color.obj and octa-color.mtl.


Many computer graphics software, however, limit the number of colors per object to a relatively small number, as 16 or 32. Blender, however, allows 32767 colors per object, a comfortable space for visualization of complex images, such as brain maps.

Stanford Polygon (*.ply)

The Stanford Polygon file allows not only the geometry to be stored, but also allow colours, transparency, confidence intervals and other properties to be assigned at each vertex. Again, the octahedron would have a file structure as:


This format also allows color attributes to be assigned to each vertex, as shown below.


The resulting object would then look like:

Batch conversion

Conversion between FreeSurfer formats can be done using mris_convert. A script to batch convert all FreeSurfer surface and curvature files from binary format to ascii is available here: subj2ascii. Before running, make sure FreeSurfer is correctly configured, that the variable SUBJECTS_DIR has been correctly set, and that you have writing permissions to this directory. Execute the script with no arguments to get usage information. The outputs (*.srf and *.dpv) are saved in a directory called ascii inside the directory of the subject. For many subjects, use a for-loop in the shell.

For more information

  • A description of the FreeSurfer file formats is available here.
  • A description of the different vtk (legacy and current) is available here.
  • Specification of the Wavefront obj format is here. Material library files are described here.
  • A description of the Stanford ply format is here.

Update — 28.May.2012: A script to convert from FreeSurfer .asc/.srf to .obj is here.
Update — 26.Jul.2013: This post was just updated to include examples of color formats.

Converting OASIS brains to NIFTI

The OASIS dataset consists of a number of T1-weighted mri brain scans, which has been kindly offered online at http://www.oasis-brains.org. The dataset can be downloaded for free after accepting an academic agreement available at the website. The data was released in analyze 7.5 file format. This format, despite having been used by the neuroimaging community for many years, suffers from not including orientation information. The nifti format addresses this concern, but when converting, some care has to be taken to ensure that there are no left-right flips. Fortunately, for this dataset, during acquisition a capsule of vitamin E was placed next to the head of each subject, on the left side, which is very helpful to identify the correct side (Marcus et al., 2007).

There are still problems, though. The nifti format has not been fully implemented in all common software packages and worse, some packages interpret differently the information contained in the header. Images that look fine in fsl‘s FSLview or FreeSurfer‘s Freeview may look stretched or shrunken in spm for instance. And images that look fine in both, may still be oriented incorrectly in Mango. Although a bit cumbersome, the procedure below ensures that the oasis images can be converted from analyze to nifti in a way that it can correctly read and shown by Mango, fsl, FreeSurfer and spm. The procedure uses exclusively fsl tools:

  1. Convert to nifti:
    fslchfiletype NIFTI_GZ OAS1_0001_MR1_mpr-1_anon.hdr
  2. Make sure there is no undesired orientation information:
    fslorient -deleteorient OAS1_0001_MR1_mpr-1_anon.nii.gz
  3. Set the sform_code as 2, which is for “aligned anatomy”. Although this is still in native, not aligned space, it ensures that software will read them appropriately:
    fslorient -setsformcode 2 OAS1_0001_MR1_mpr-1_anon.nii.gz
  4. Set the sform as the following matrix:
    fslorient -setsform  0 0 -1.25 0  1 0 0 0  0 1 0 0  0 0 0 1  OAS1_0001_MR1_mpr-1_anon.nii.gz
  5. Swap the order of the data. Again, this isn’t really necessary, except to ensure that different applications will all read correctly:
    fslswapdim OAS1_0001_MR1_mpr-1_anon.nii.gz RL PA IS OAS1_0001_MR1_mpr-1_anon.nii.gz
  6. fsl tries to preserve orientation and, when the voxels are reordered, it modifies the header accordingly, resulting in no net transformation when seen with fsl tools. To resolve this, it’s necessary to change the header again, now the qform:
    fslorient -setqform -1.25 0 0 0  0 1 0 0  0 0 1 0  0 0 0 1  OAS1_0001_MR1_mpr-1_anon.nii.gz

These steps can be placed inside a simple loop within the shell for either Linux and Mac, like below (click here to download):

#!/bin/bash

# Directory where the OASIS data is located
ROOTDIR=/Volumes/HD2/oasis-dataset-416subj/original

# Directory to save results
NIFTIDIR=${ROOTDIR}/../nifti

# Go to the directory with the data
cd ${ROOTDIR}

# For each subject
for s in * ; do

  # Some feedback in the screen
  echo ${s}

  # Create directory to save the results, if not existing
  mkdir -p ${NIFTIDIR}/${s}

  # Directory of the original, raw data
  cd ${ROOTDIR}/${s}/RAW

  # For each acquisition
  for a in *.hdr ; do

    # Do each of the 6 steps described in the blog
    ${FSLDIR}/bin/fslchfiletype NIFTI_GZ ${a} ${NIFTIDIR}/${s}/${a%.hdr}
    ${FSLDIR}/bin/fslorient -deleteorient ${NIFTIDIR}/${s}/${a%.hdr}
    ${FSLDIR}/bin/fslorient -setsformcode 2 ${NIFTIDIR}/${s}/${a%.hdr}
    ${FSLDIR}/bin/fslorient -setsform  0 0 -1.25 0  1 0 0 0  0 1 0 0  0 0 0 1  ${NIFTIDIR}/${s}/${a%.hdr}
    ${FSLDIR}/bin/fslswapdim ${NIFTIDIR}/${s}/${a%.hdr} RL PA IS ${NIFTIDIR}/${s}/${a%.hdr}
    ${FSLDIR}/bin/fslorient -setqform -1.25 0 0 0  0 1 0 0  0 0 1 0  0 0 0 1  ${NIFTIDIR}/${s}/${a%.hdr}
  done
done

echo "Done!"

The reference for the oasis dataset is: