# Retrospective motion correction for structural MRI

Good quality results can only be obtained with good quality images. To obtain images with good SNR and high resolution, long scan times are often necessary. However, with longer acquisition, the subjects are more likely to move. Rather than acquire a very long sequence, that often can last up to 40 minutes, the acquisition can be broken into many shorter scans, each with lower SNR. These are coregistered, effectively reducing the effect of undesired movement, and averaged, for increase of the SNR (Kochunov et al, 2006).

The last image in the second row is the average of all the other seven, after retrospective motion correction has been applied.

The results above are typical. Each pixel above corresponds to a true voxel, with no interpolation. Therefore, to see in more detail, simply download the image above and use an image viewer (e.g. gwenview or eog) to zoom.

The method is termed retrospective to discern it from prospective motion correction methods, such as those available for BOLD-sensitive EPI acquisitions (e.g. Siemens 3D PACE).

Since the images belong to the same subject and typically the movement is of only a few degrees or milimeters, the affine matrices for these corrections should have a diagonal with values very close to 1 and very close to zero off the diagonal. Here is one example:

0.999995    -0.00020007  -0.0030414    0.106906
0.000188555  0.999993    -0.00378601   0.100499
0.00304214   0.00378542   0.999988    -0.204133
0            0            0            1

As usual, the improvement in the SNR is given by a factor of $\sqrt{n}$, where $n$ is the number of images averaged.

A script to apply the method is available here: stereomc. Just download and make it executable. Call it without arguments to obtain usage information.

## Requirements

FSL tools must be installed and properly configured.

## Reference

Kochunov P, Lancaster JL, Glahn DC, Purdy D, Laird AR, Gao F, Fox P. Retrospective motion correction protocol for high-resolution anatomical MRI. Hum Brain Mapp. 2006 Dec;27(12):957-62.

# The logic of the Fisher method to combine P-values

Consider a set of $k$ independent tests, each of these to test a certain null hypothesis $\mathcal{H}_{0|i}$, $i=\{1, 2, \ldots, k\}$. For each test, a significance level $P_{i}$, i.e., a p-value, is obtained. All these p-values can be combined into a joint test whether there is a global effect, i.e., if a global null hypothesis $\mathcal{H}_0$ can be rejected.

There are a number of ways to combine these independent, partial tests. The Fisher method is one of these, and is perhaps the most famous and most widely used. The test was presented in Fisher’s now classical book, Statistical Methods for Research Workers, and was described rather succinctly:

When a number of quite independent tests of significance have been made, it sometimes happens that although few or none can be claimed individually as significant, yet the aggregate gives an impression that the probabilities are on the whole lower than would often have been obtained by chance. It is sometimes desired, taking account only of these probabilities, and not of the detailed composition of the data from which they are derived, which may be of very different kinds, to obtain a single test of the significance of the aggregate, based on the product of the probabilities individually observed.
The circumstance that the sum of a number of values of $\chi^{2}$ is itself distributed in the $\chi^{2}$ distribution with the appropriate number of degrees of freedom, may be made the basis of such a test. For in the particular case when $n=2$, the natural logarithm of the probability is equal to $\frac{1}{2}\chi^{2}$. If therefore we take the natural logarithm of a probability, change its sign and double it, we have the equivalent value of $\chi^{2}$ for 2 degrees of freedom. Any number of such values may be added together, to give a composite test, using the Table of $\chi^{2}$ to examine the significance of the result. — Fisher, 1932.

The test is based on the fact that the probability of rejecting the global null hypothesis is related to intersection of the probabilities of each individual test, $\prod_{i} P_{i}$. However, $\prod_{i} P_{i}$ is not uniformly distributed, even if the null is true for all partial tests, and cannot be used itself as the joint significance level for the global test. To remediate this fact, some interesting properties and relationships among distributions of random variables were exploited by Fisher and embodied in the succinct excerpt above. These properties are discussed below.

## The logarithm of uniform is exponential

The cumulative distribution function (cdf) of an exponential distribution is:

$F(x)=1- e^{-\lambda x}$

where $\lambda$ is the rate parameter, the only parameter of this distribution. The inverse cdf is, therefore, given by:

$x = -\dfrac{1}{\lambda}\ln(1-F(x))$

If $P$ is a random variable uniformly distributed in the interval $[0, 1]$, so is $1-P$, and it is immaterial to differ between them. As a consequence, the previous equation can be equivalently written as:

$x = -\dfrac{1}{\lambda}\ln(P)$

where $P \sim \mathcal{U}(0,1)$, which highlights the fact that the negative of the natural logarithm of a random variable distributed uniformly between 0 and 1 follows an exponential distribution with rate parameter $\lambda=1$.

## An exponential with rate 1/2 is chi-squared

The cdf of a chi-squared distribution with $\nu$ degrees of freedom, i.e. $\chi^{2}_{\nu}$, is given by:

$F(x; \nu) = \dfrac{\int_{0}^{x/2} t^{\frac{\nu}{2}-1}e^{-t}{\rm d}t}{\left(\frac{\nu}{2}-1\right)!}$

If $\nu=2$, and solving the integral we have:

$F(x; \nu=2) = \dfrac{\int_{0}^{x/2} t^{\frac{2}{2}-1}e^{-t}{\rm d}t}{\left(\frac{2}{2}-1\right)!} = \int_{0}^{x/2} e^{-t}{\rm d}t = 1-e^{-x/2}$

In other words, a $\chi^{2}$ distribution with $\nu=2$ is equivalent to an exponential distribution with rate parameter $\lambda=1/2$.

## The sum of chi-squared is also chi-squared

The moment-generating function (mgf) of a sum of independent variables is the product of the mgfs of the respective variables. The mgf of a $\chi^{2}_{\nu}$ is:

$M(t) = (1-2t)^{-\nu/2}$

The mgf of the sum of $k$ independent variables that follow each a $\chi^{2}_{2}$ distribution is then given by:

$M_{\text{sum}}(t) = \prod_{i=1}^{k} (1-2t)^{-2/2} = (1-2t)^{-k}$

which also defines a $\chi^{2}$ distribution, however with degrees of freedom $\nu=2k$.

## Assembling the pieces

With these facts in mind, how to transform the product $\prod_{i} P_{i}$ into a p-value that is uniformly distributed when the global null is true? The product can be converted into a sum by taking the logarithm. And as shown above, the logarithm of uniformly distributed variables follows an exponential distribution with rate parameter $\lambda=1$. Multiplication of each $\ln(P_{i})$ by 2 changes the rate parameter to $\lambda=1/2$ and makes this distribution equivalent to a $\chi^{2}$ distribution with degrees of freedom $\nu=2$. The sum of $k$ of these logarithms also follow a $\chi^2$ distribution, now with $\nu=2k$ degrees of freedom, i.e., $\chi^{2}_{2k}$.

The statistic for the Fisher method is, therefore, computed as:

$X = -2 \sum_{i=1}^{k} \ln(P_{i})$

with $X$ following a $\chi^{2}_{2k}$ distribution, from which a p-value for the global hypothesis can be easily obtained.

## Reference

The details above are not in the book, presumably omitted by Fisher as the knowledge of these derivation details would be of little practical use. Nonetheless, the reference for the book is:

The Fisher’s method to combine p-values is one of the most powerful combining functions that can be used for Non-Parametric Combination.

# Importing FreeSurfer subcortical structures into Blender

In a previous article, a method to import cortical meshes from FreeSurfer into Blender and other CG applications was presented. Here it is complemented with a method to import the subcortical structures. These structures are segmented as part of FreeSurfer subcortical stream, and are stored in a volume-based representation. In order to import into Blender, first it is necessary to generate surfaces, then convert to Wavefront OBJ, and finally, do the import. Below is an example of how the output looks like.

An example of a subcortical meshes generated from FreeSurfer subcortical segmentations imported into Blender.

## Generating the surfaces for subcortical structures

The surfaces can be generated using FreeSurfer’s own tools. The script aseg2srf automates the task. Suppose you would like to create meshes for all subcortical structures of the subjects named “capitu”, “bentinho” and “ezequiel”. Simply run:

aseg2srf -s "capitu bentinho ezequiel"


# Importing FreeSurfer cortical meshes into Blender

Computer graphics (CG) applications, such as Blender, Maya or 3D Studio, doesn’t come with filters that allow to directly import or export FreeSurfer meshes. To import brain surfaces to these applications, converting first to an intermediate file format solves the problem. Below is an example of a brain imported into Blender.

An example of a brain mesh from FreeSurfer imported into Blender after conversion to Wavefront OBJ.

The mris_convert, from FreeSurfer, allows conversion to stereolithography format (STL). However, not all CG applications read natively from STL files. Moreover, other formats, such as Wavefront OBJ and Stanford PLY have simpler structure — yet allowing features that will be discussed here in forthcoming articles.

## Conversion from FreeSurfer to Wavefront OBJ

A simple procedure to convert a FreeSurfer surface to Wavefront OBJ (ASCII) is:

1) Use mris_convert to convert from FreeSurfer binary to ASCII:

mris_convert lh.pial lh.pial.asc


2) A suggestion is to rename the .asc file to .srf, to differentiate it from other files that also have .asc extension and have a different internal structure.

mv lh.pial.asc lh.pial.srf


3) Use the script srf2obj to convert from .srf to .asc:

srf2obj lh.pial.srf > lh.pial.obj


The symbol “>” is necessary so that the output from srf2obj is redirected to a new file.

## Conversion from Wavefront OBJ to FreeSurfer

To convert back to FreeSurfer, a similar strategy is used.

1) Use the script obj2srf to generate new .srf files (here directly named as .asc, as expected by mris_convert)

obj2srf lh.pial.obj > lh.pial.asc


2) Use mris_convert:

mris_convert lh.pial.asc lh.pial