ICA-based artifact removal diminishes scan site differences in multi-center resting-state fMRI

Resting-state fMRI (R-fMRI) has shown considerable promise in providing potential biomarkers for diagnosis, prognosis and drug response across a range of diseases. Incorporating R-fMRI into multi-center studies is becoming increasingly popular, imposing technical challenges on data acquisition and analysis, as fMRI data is particularly sensitive to structured noise resulting from hardware, software, and environmental differences. Here, we investigated whether a novel clean up tool for structured noise was capable of reducing center-related R-fMRI differences between healthy subjects. We analyzed three Tesla R-fMRI data from 72 subjects, half of whom were scanned with eyes closed in a Philips Achieva system in The Netherlands, and half of whom were scanned with eyes open in a Siemens Trio system in the UK. After pre-statistical processing and individual Independent Component Analysis (ICA), FMRIB's ICA-based X-noiseifier (FIX) was used to remove noise components from the data. GICA and dual regression were run and non-parametric statistics were used to compare spatial maps between groups before and after applying FIX. Large significant differences were found in all resting-state networks between study sites before using FIX, most of which were reduced to non-significant after applying FIX. The between-center difference in the medial/primary visual network, presumably reflecting a between-center difference in protocol, remained statistically significant. FIX helps facilitate multi-center R-fMRI research by diminishing structured noise from R-fMRI data. In doing so, it improves combination of existing data from different centers in new settings and comparison of rare diseases and risk genes for which adequate sample size remains a challenge.

section_content: Registering scans of brains from different subjects is a neccesary processing step for many types of analyses such as e. g. multi-subject fMRI studies, inter-group comparisons of tissue composition (e. g. VBM) or measures derived from diffusion weighted MR (e. g. TBSS). There is a family of algorithms that attempt to do this in the native "brain-space" using different sets of transforms between subjects, or more typically between a subject and some template instantiating some standard space. In order of allowing more "non-local" warps the transforms can be divided into linear (affine) transforms, small-deformation non-linear transforms and large-deformation non-linear transforms (often referred to as "viscous fluid registration"). The principal difference between small-deformation and large-deformation method is that in the former case the warps are fully described by three fields of displacements (one for each dimension). In the fluid-type algorithms one would in principle need the history of each location as it is being warped towards it target registration, though in practice this is often achieved through one or more regriddings within the registration. An intuitive example of the difference between the approaches is the matching of two 2D pictures of human faces. For small-deformation methods it would e. g. not be possible to map the left eye of subject one to the right eye of subject two and vice versa. This is in contrast to large-deformation methods where it would in principle be possible to start by displacing the right eye upwards and the left downwards, move each eye towards the opposite side and then move them up/down into place. 

The fluid-registration algoritms are often considered the most "advanced" and are typically formulated in terms of a set of partial differential equations, and are often solved using corresponding methods such as Gauss-Seidel or multi-grid methods. As outlined above they offer great freedom in matching one image to another and can achieve an almost arbitrarily good matching of intensities between the images. However, the information available in "traditional" structural images such as T1-or T2-eighted images is largely limited to tissue type. This means that prety much any mapping between two images that map gray matter onto gray matter, white matter onto white matter and CSF onto CSF is equally good in terms of the cost-function and the task of finding the "true" warps can be thought of as finding the most likely set of warps of those that yield that best cost-function value. It can therefore be argued that warps from a small-deformation model that yields a similar or identical cost-function value as those from a large-deformation model represents a more "reasonable" solution to the registration problem. Furthermore, it is clear that there exist no one-to-one mapping as defined by gyri and sulci between different brains. There is a set of sulci that are consistently found across "normal" subjects such as e. g. the interhemispheric fissure, the Sylvian fissure, the parietal-occipital fissure and the central sulcus. In contrast there appear to be large inter-subject differences in e. g. parietal cortex where it has been reported that not even the number of sulci is consistent across "normal" subjects. Even though it would seem reasonable to assume that there exist at least a functional one-to-one mapping across subjects there is insufficient experimental data to conclude even that. Indeed there exists studies which indicate that even low-level tasks such as e. g. odd-ball detection is processed differently in different subjects. 

In this work we have therefore opted for a small-deformation model for the warps, and put our efforts into ensuring convergence to a plausible field. We have aimed at achieveing that through a combination of a pyramid, or multi-resolution, scheme along with an afficient optimisation algorithm for each step of the pyramid.
section_name: INTRODUCTION

section_name: THEORY

section_content: The general transform of coordinates for nonlinear registration is of the form

where M is an affine transform matrix that will account for differences in voxel-size, position, etc. For the reminder of this paper we will assume that M is known. The entities d i (x, y, z) are knowns as displacement fields, and will for each location [x y z] tell how far the sampling point should be displaced in the i direction. 

We will further assume the existence of some interpolating function such that given samples of some function g(i, j, k) where i, j and k are sampled on some regular grid we can infer values g(x, y, z) as long as x, y and z fall within the original grid. We will then use g(x, y, z ) to denote the values of g for a set of coordinates given by 1. This value can be thought of as the value of g at some point [x y z] under the transform given by M and the displacement fields d i (x, y, z). 

In this paper we will use basis-functions to model the fields d i (x, y, z) as a function of some set of parameters w. We can for example model the field as a linear combination of 3D cubic B-splines. Each spline in the 3D set is associated with a spline coordinate lmn specifying the position of the spline in the space given by xyz. We will use B lmn (x, y, z) to denote the value of the spline with coordinate lmn at the location [x y z]. Each of these splines are then multipled by a coefficient c lmn, specifying how much is needed of that particular spline (see fig. ?). We now vectorise the coefficients c lmn and denote the resulting vector by w i where i denotes the direction of the field such that e. g. w x denotes the vector containing the coefficients defining the d x (x, y, z) displacement field. Furthermore we denote the concatenation of these vectors by w, i. e w = [w (x)T w (y)T w (z)T ] T. With this notation we can the write g xyz (w) to denote the value of g at [x y z ] under the transform given by w implicit on the (constant) M. 

The general notation g xyz (w) is independent of our choice of basis-function and would be equally valid for e. g a discrete cosine basis set. 

We will also make use of the gradient, or rather its components, of g at some point [x y z ], which we define as

where the partial derivative with respect to e. g x, ∂g xyz /∂x| w, denotes the rate of change of g at [x y z ] as one translates the sampling point in the x-direction. It should be noted that the exact form of the partial derivatives ∂g/∂x is conditional on the interpolating function alluded to above. In the present paper we will ignore this issue and assume the existence of functions/kernels to perform the interpolation and calculating the corresponding partial derivatives. If we now assume that the parameter w

is the coefficient for the lmnth spline in the x-displacement field d x we can denote the derivative of g as

where the w subscript indicates that is has been calculated at a point w in the parameter space. 

Let us further define the vector

where x, y and z are defined on the ranges 1-X, 1-Y and 1-Z respectively. In an equivalent manner we define the vectors

and

Noting that there is a direct mapping from the index i into w to the triplet lmn we can equally well, and more succinctly, write B i as B lmn for the vector representation of the ith basis function. 

Combining equations 3, 5 and 6 we can further define the matrix

where denotes elementwise (or Hadamard) product. Hopefully it should be clear that each element of J x (w) is of the form given by equation 3. Finally by thinking of g(w) as an 3LM N → XY Z mapping where L, M and N are the number of basis-functions and X, Y and Z are the number of samples/voxels in the x-, y-and z-directions respectively we can define the Jacobian matrix J of that mapping as
section_name: Transforms and some definitions

section_content: Our task is to find the parameters w describing the fields d i (x, y, z) that transforms our function/image g from its native space to some other arbitrary reference space. In this paper we consider intensity-based methods for finding these. This is performed by defining some cost/objective-function in terms of some template f that defines our reference space and our object image g(w). Let us call this function O, and let us write it as O(w) as a short for O(f, g(w)). 

An example of such a function is the "mean sum of squared differences" cost-function which we can write as

or equivalently using the notation we defined in equation 4

where the subscript w indicates that the derivative has been calculated at that point in the parameter space. This is the transpose of how it is mostly defined in the literature, but it will save us from writing some T s in the remainder of the paper. From the definition of O in equation 9 we see that an element ∂O/∂w j of ∇O can be written as

If we define a vector e as

equivalently to equations 4-6, and using the definition in equation 8 we see that the gradient of O can be written as

The next entity that we are intersted in is the Hessian of O, i. e. the matrix whose ijth element is

Again from the definition of O given in equation 9 we see that such an element is of the form

Using equations 8 and 17 we see that we can write the Hessian matrix H as

J T (w)J(w) + term with weird second derivatives (18)

The first term in equation 18 is known as the Gauss-Newton approximation to the Hessian. In our technical report on optimisation we explain why this works, and often works better than if using the exact Hessian. For the remainder of this paper we will work with an approximation of the hessian given by
section_name: Sum of squared differences cost-function

section_content: In the implementation of fnirt there is a virtual base class called basisfield from which two classes, splinefield and dctfield, have been derived. An object of such a class contains information about the field it implements. So if we were to e. g to declare a field of type splinefield we would do it like Thus we have an object called xfield, corresponding to the entity d x (x, y, z) in the text above, implementing a cubic B-spline field of size 128 × 128 × 96, corresponding to X × Y × Z in the equations above, with a knot-spacing of 8 voxels. The latter means that the splines are placed on a regular grid with a distance of 8 voxels between adjacent spline kernels. These are placed such that the centre of spline # 2 coincides with the centre of the first voxel and then placed at regular intervals of 8 voxels until the first kernel whose whole support falls outside the field. From this are given the entities L, M and N in the equations above, and we can enquire about these of our field as Let us now say we have vectors f, g, and dgdx containing the reference image f, the object image g and the partial derivative w. r. t. x, ∂g/∂x, both sampled at some point w in reference space. And let us say the we now want to calculate the gradient and Hessian of O. We can then note that We want of course all of nablaO and we would get that from corresponding field objects for d y (x, y, z) and d z (x, y, In the code above we could actually make do with a single field, provided we want to model the x-, y-and z-displacement fields with the same resolution, which we assume we always want. The important thing in the example above is that the member function Jte gets called with different partial derivatives. However, in any actual case we will still need to represent the three displacement fields. 

The next entity we would like to calculate is the (approximate) Hessian of O. Looking at equations 8 and 19 we see that it can be written as

Let us now look at how we would calculate two of the submatrices, J T x (w)J x (w) and J T x (w)J y (w). Please don't be confused by the declaration of Hxx and Hxy. It is basically a pointer to an object of type BFMatrix. The class BFMatrix in turn is a virual base class with two derived classes SparseBFMatrix and FullBFMatrix which implements interfaces to sparse and full matrices respectively. 

The final thing that should be pointed out in this section is that it is very easy to write code that is independent of what basis-set we actually want to use. The following code would for example work So, you can see that after the initial creation of the field we can do all we need to do with it without actually having to know what type of field it is.
section_name: Relation to my code

section_content: This is a tiny change compared to the "sum of squared differences" described above. But it is of interest because it is the first (most basic) actual cost-function implemented in fnirt, and as such will allow us further insights into the structure of the code. For the more general cases later on it is now useful to define the vector θ as the concatenation of w, the parameters pertaining to the displacement field, and any other parameters we might want to introduce. These other parameters will typically be used to model variations/differences in intensity that does not pertain to structure and can hence be seen as confounds within the present problem. The very simplest form would be a scalar scaling factor that models global intensity differences between f and g. So, let us now define θ as

and the cost-function as

which yields the gradient

and the (approximate) Hessian

3. 4. 1 So, how do we calculate that using the basisfield class?

The upper left block is identical to before, so we have already seen how to calculate that. It can further be seen that the block consisting of -J T (w)f is of the same general form as the upper portion of the gradient, only with the vector (g(w) -αf ) replaced by f. As tempting as it might seem to try this out straight away, we should mention that there is yet another layer of abstraction though.
section_name: Sum of squared differences cost-function with scaling

section_content: As part of implementing a toolbox for nonlinear optimisation we defined a cost-function class with the very simple interface It can be seen that the member function cf is pure virtual, which means that NonlinCF is a virtual base class. Its purpose is to provide a consistent interface for obtaining the cost-function, its gradient and Hessian. The parameter p corresponds to θ in the paragraph above and e. g a call like mygrad = cf. grad(p) is literally identical to ∇O(theta). 

Being a virtual base class there will never be any instances of NonlinCF but rather of classes derived from it. To make it really concrete let us imagine we want to fit a mono-exponential function to some data under the assumption that errors are normal distributed. So, our model is given by

where we are interested in finding θ = [θ 1 θ 2 ] T. We do so by defining a cost-function which is the sum of squared errors between the model predictions and the observed data. As you see we have now added "space" for the data in the class, and defined the member function cf. We are a little lazy so we don't override neither grad or hess. This is OK since the base class NonlinCF defines these using numerical differentiation based on the cf function that we have just defined. To use this we will further need to create an instance of the class NonlinParam which is a glorified struct that contains information about what algorithm we want to use for the minimisation, convergence critera etc. So, all the code we need to write is As can be deduced from the above nonlin is a global function that takes as parameters a NonlinPar and a NonlinCF (or any class derived from NonlinCF) object and finds a set of parameters that minimise the value of the cost-function. There is a rich set of parameters that can be set for objects of type NonlinParam that should be sufficient to ensure convergence for most types of models and data.
section_name: The cost-function class and the nonlinear toolbox

section_content: The fnirt CF class is intended as a base-class from which to derive cost-function classes for use in fnirt (or fnirt-like applications). It is derived from NonlinCF, but contrary to what its name implies it doesn't actually implement any cost-function as such. Its purpose is instead to implement functionality that transcends any particular cost-function. It does so through a set of protected member functions (i. e. functions that are only available from within subclasses of fnirt CF). It may all seem a little esoteric so let us dive straight into an example So, making ones cost-function a subclass of fnirt CF makes it relatively easy to implement the cost-function. The functionality inherited from NonlinCF makes it easy to perform the actual minimisation and the functionality from fnirt CF facilitates implementing the actual cost-function (and its derivative and Hessian). So in the example above we used a mask in object space because we didn't want the warping to be confused by some ghastly growth that happens to be in poor obj. Because the mask is in the space of the image we attempt to warp it means that it will have to be warped along with obj in each iteration. This is facilitated by fnirt CF offering a public interface. SetObjMask(mask) that allows us to specify the mask and the protected function. Mask() that returns a mask that is the intersection of all the masks (in object and/or reference space) that has been set by the user transformed into reference space. This is just one example of how the functionality in fnirt CF may facilitate the implementation of new cost-functions.
section_name: The fnirt CF class

section_content: The chosen model for non-linear registration offers considerable freedom in terms of different constellations of warps. It is frequently the case that a large set of different warps (as parametrised by different instances of w) yields similar values for the cost-function, and in these cases we need some way to decide between them. Furthermore, representing the field as a linear combination of basis-fuctions ensures that it is smooth and continous but does not not guarantee that it is "one-to-one" and "onto". 

The term "one-to-one" means that no two points in the original space x can map onto the same point in the warped space x, a condition which is indicated by the Jacobian determinant of the [x y z] → [x y z ] transformation becoming zero or negative at one or more grid points. The significance of "onto" is that each point in the transformed space x should have a mapping onto some point in the original space x, i. e. there must not be any points in the transformed space that "cannot be reached". The "onto" requirement is not really meaningful (i. e. we can never hope to enforce it) for a discretely sampled function, but the "one-to-one" condition is widely thought to be important if we are to consider a transform/field as reasonable. 

As alluded to above the "one-to-one" requirement is not guaranteed, or even helped, by representing the field by some basis-set. All that guarantees is that when the Jacobian goes negative it does so gradually (over space). A common solution to these problems is to use some form of "regularisation" on the field. This is simply some differentiable function of the field, or the parameters of the field, whose value indicates how "likely" we consider that field to be. It is quite common to use some mechanical analogy such that in the choice between two fields we will consider e. g. that with a smaller membrane energy to be the more likely. The mechanical analogy results in the same set of equations as one obtains from considering the sum-of-squared differences cost-function as a likelihood (assuming normal distributed errors) and postulating a multinormal prior on the coefficients of the field. There are forms for the variance-covariance matrix of the prior that corresponds to e. g. membrane energy and bending energy. 

At present the membrane energy is used for regularisation. Given the definition in equation 1 the membrane energy is defined as

where the i subscript denotes the ith voxel and where λ is a material dependent constant. If we define a vector B (i) lnm as a vector containing an "unravelled" version of the lmnth basisfunction having been differentiated in the ith direction (cf equation 6), the matrix B (i) as the concatenation of the LM N basis-functions (cf equation 7) and the vector w x containing the coefficients pertaining to the x-component (d x (x, y, z)) of the field then the membrane energy for the x-component can be expressed as

and the total membrane energy as

Consequently the gradient and the Hessian of E m (w) are

and

Hence the cost-function including regularisation, its derivatives and Hessian are

,

and

where λ can be seen as an arbitrary weighting of the regularisation versus the sum-of-squared differences, or between the likelihood and the prior if one prefers that perspective.
section_name: Regularisation of the field

section_content: The contribution of the membrane energy is a function only of the field, or rather the parameters defining the field, and can hence be calculated and returned by the basisfield class. Let us for example say that we have declared and defined a splinefield like
section_name: Implementation in the basisfield class

section_content: The Hessian H will rapidly become very large and costly to store and calculate (more about the latter below). A reasonable "bread-and-butter" resolution may be, in the case of using a spline-basis, a 4 × 4 × 4 voxel knot-spacing over the 91 × 109 × 91 matrix of the MNI template. This would result in 25 × 30 × 25 splines requiring 18750 parameters for each of the x-, y-and z-displacement fields. This means that H is a 56251 × 56251 matrix, requiring on the order of 25GB to store and represent it, something which is beyond most computers today. This is the reason why for example the DCT-based registration in SPM is limited to a resolution of roughly 20mm (or 10 voxels) isotropically. The spline basis set has a great advantage here in that the Hessian is sparse, meaning that a large proportion of the entries are zero. And that proportion increases with increased resolution of the warps. No row or column of the Hessian contains more than 343 (7 × 7 × 7) non-zero elements, and many contain fewer. It is therefore crucial that we have an efficient storage format for the Hessian (or for sparse matrices in general) if we are to reap the full benefits of the spline basis-set. At the same time we do not wish to represent all Hessian matrices as sparse matrices (e. g. the full Hessian matrices from the DCT basis set) since that will typically increase storage needs by more than 50% and incurr an execution time penalty when the matrix is indeed full. This is the reason why our basisfield classes, and also the classes derived from NonlinCF, return an object of type BFMatrix. This is a virtual wrapper-class with two derived classes FullBFMatrix and SparseBFMatrix has an API that is sufficient for the needs of the different sub-classes of NonlinCF. This is a way to obtain the required polymorphism given the restriction that the matrix-class we commonly use, Newmat, do not have a sparse representation nor has it been designed such that it would be easy to sub-class a sparse-matrix class from it. 

Why do we then bother to include the DCT basis-set in the first place? When performing low-resolution ( 20mm warp resolution or poorer) non-linear warping the full Hessian for the DCT set is feasible to represent and very efficient to calculate. This has to do with the properties of the basis and slightly simplified it can be though of as an FFT for Hessian calulation. In contrast the Hessian for the spline set is not very sparse at that resolution, and each element is quite costly to calculate because of the large overlap (in units of voxels) of neighbouring splines. The DCT-set will therefore yield similar results as the spline-set for an order of magnitude shorter execution time. 

The advantage of the spline-set becomes obvious when going to medium-resolution registration ( 10mm warp resolution), when it is no longer possible to represent the full Hessian for the DCT-set. Even if it was possible to represent it the Hessian contains 64 times more elements, each element taking as long to compute as for the low-resolution case. For the spline set the number of elements grow much more slowly with increasing reolution since sparsity also increases with resolution. In addition the cost of calculating each element decreases since the overlap between neighbouring splines decreases. In fact there is but limited difference in the time it takes to calculate the 5616 × 5616 Hessian for a knot-spacing of 20mm and the 31752 × 31752 Hessian for a 10mm knot-spacing.
section_name: What is going on with that BFMatrix class?

section_content: When using Gauss-Newton style optimisation (Levenberg or Levenberg-Marquardt) the costliest operations are the calculation and inversion of H. The example in the preceeding section shows that a 56251 × 56251 Hessian is in no way unreasonable, and will take a little while both to compute and invert. For this reason the. JtJ() routines of the derived classes of basisfield are the real bottle necks, and efforts at optimisation should focus on these. 

If we look at e. g. equation 19 we see that H is a function of w and will thus need to be recalculated at each iteration. If we further look at equations 7 and 8 we see that that dependence comes from caculating ∂g/∂x (and y and z) at the point w in the parameter space. What then do we think ∂g/∂x would look like when we are reasonable close to convergence? I think it would look very similar to ∂f /∂x, and that doesn't depend on w. So one trick, that has been played successfully in e. g. SPM, is to calculate H once and for all using ∂f /∂x. The approximation is of course worse for the first few iterations when g and f are still quite dissimilar. However, at that stage the second order Taylor expansion that underlies Gauss-Newton style minimisation is likely to be a poor approximation anyway. It should also be noted that getting H "wrong" in a Levenberg/Levenberg-Marquardt minimisation should not affect the end result, just the way and time it takes to get there. A poor approximation to H would mean that would take more iterations to converge, but on the other hand that should be weighed against the time and computational effort of each iteration which will be much smaller if we can reuse H. Even better, of course, would be if we could in the first iteration calculate the Cholesky decomposition of H which would greatly facilitate solving for θ in subsequent iterations. 

These are all things that we need to look into empirically. It is neither clear that we will end up using Gauss-Newton style optimisation as our default. In particular the Scaled Conjugate Gradient method seems to do well on the model examples I have looked at, and that does away with calulating H altogether (though it calculates the gradient twice at each iteration). However, see my reservations below regarding different scales of different sets of parameters.
section_name: Approximations for speed

section_content: Many algorithms for non-linear optimisation suffers quite badly when different parts of the gradient vector (∇O) have vastly different scales. Consider ∇O given by equation 25. It is a 3LM N + 1 × 1 vector where the upper 3LM N elements are of the form j T i (g(w) -αf ) where j i is a vector which is non-zero only over small portion corresponding to the support of spline that it pertains to. This is in contrast to the final element -f T (g(w) -αf ) where more or less all elements of f are non-zero. Intuitively it is also obvious that a unity change of α (say from 1 to 2) will change O massively, whereas a unity change of one of the spline coefficients will move the samling point for a few hundred voxels to varying degrees and will have much smaller effect on O. So we have a vector where 3LM N elements are of one scale, and then one element that is hundreds, or even thousands, of times larger. 

Let us then consider what would happen when a method that bases its step (in parameter space) on the gradient, i. e. attempts to take a step in the -∇O(θ) direction. That direction will be completely dominated by the α direction, i. e. we attempt to take a long step in the α direction and a short step in all other directions. The problem though is that most likely we only need to take a very short step in the α direction. If we have done an intial global normalisation (as one does) a realistic range might be 0. 97 < α < 1. 03. So that means that for the step θ (k+1) = θ (k) -λ∇O(θ) we will find a very small λ, and therefore it is likely that we will only have taken a tiny step of that we would need to take for the other direction. We should note here that although the search directions for e. g. Variable-Metric and Conjugate-Gradient methods are not the gradient, they both start ot with that as their initial direction. And because of the way they construct new directions based on previous directions and the present gradient each new directions will have a large proportion of the present gradient in them. 

In practice what will tend to happen is that for the first few (can be a small few, but sometimes a large few) attempted steps the gradient will be dominated by -f T (g(w) -αf ) so the step lengths become very small and we practicaly do not move at all in any of the w directions. Then comes a step where we are literally at the bottom of the pass in the α direction and ∂O/∂α becomes (very close to) zero. The algorithm will then be able to take a sizeable step along the w directions and the cost-function will decrease. After that step α is typically non-zero again and there will be another few steps of neglible length until ∂O/∂α is again zero and the algorithm will again be able to take a "proper" step. So, the problem produces very typical symptoms of slow and jerky convergence. For the particular case of a single scale-factor whose scale is vastly different from the warp (w) it is relatively easy to "fudge" it by redefining O as a summation over (g xyz (w) -(1 + cα)f xyz ) 2 where c is some suitably small constant (10 -3 seems to work reasonably well). However the problem is more "general" than that and for example for my top-down-bottm-up method the parameter vector θ = [w T p T ] T where w is a set of spline-coefficients used to model the (magnetic) field and p are a set of rigid body movement parameters describing any change in subject position between the two acquisitions. The effects on the cost-function is of course much larger for any of the movement parameters than for the warp parameters, and we have again a problem of different scales of parameters. We will see similar problems also if/when we want to include a model for RF inhomogeneity. I suspect/hope it shall be possible to find reasonable fudge factors also for those parameters, but it won't be pretty. It is hardly realistic that we shall be able to find a general solution to it since this is a long standing issue in "non-Newton" non-linear optimisation. 

In contrast, for the Newton/Gauss-Newton style algorithms this is a non-issue. If we take again the example of the gradient described by equation 25, remember that the step for a Newton style algorithm is of the general form θ (k+1) = H(θ) -1 ∇O(θ) and look at the element in the bottom-right corner of H we see that that, just as the large element of ∇O, that too is much larger than the other elements of H. A, very crude, approximation to the α component of the Newton-step is -f T (g(w) -αf )/f T f, and in practice this typically turns out to be one of the smallest components of the step. The information in the Hessian will for all practical purposes "calibrate" any akward scalings in the gradients without any need for finding ones own fudgy factors. This means that convergence tend to be steady and reliable, which makes it easy to implement convergence criteria. Non-Newton type algorithms, as described above, will for poorly scaled parameters exhibit long stretches of sometimes tens of iterations during which changes in cost-function are very small, and may then all of a sudden make another "jerk" and change by maybe a factor of two (I have seen it happen). How, then, does one implement a test for convergence?

Intrestingly the Variable-Metric methods are actually supposed to be able to handle this type of problem (hence the name), but in my practical tests I have found nothing to indicate that would actually be the case. 

For this reason I tend to favor Gauss-Newton type algorithms over the various flavours of Variable-Metric/Conjugate-Gradient methods out there that doesn't require calculation or representation of the Hessian. Although the price of course is that it complicates the implementation quite considerably and represents a large proportion of the total calculations.
section_name: Scaling of the gradient vector

section_content: Experiments have been performed registering different structural images against each other and against templates acquired with similar tissue contrasts. This has been performed both for T1-weighted and FA (Fractional anisotropy) images. Tests have been aimed at finding Figure 1 : Sagittal slice at y=6mm in the MNI-space. On the left is an average across 10 subjects after affine and on the right after non-linear registration. The parietal-occipital fissure was manually traced in each subjects normalised scan and all ten traces were overlayed on the average images. Not the considerably higher level of detail in the non-linear mean, indicating better alignment of structure, and the considerably smaller inter-subject variability of the manual traces. 

schemes/schedules of combinations of subsampling of images and field, regularisation weight and smoothing of images.
section_name: EXPERIMENTS

section_content: An example of the results from the non-linear registration is demonstrated in figures 1 and 2. These are T1-weighted structural images that were registered to the non-linear avg152 template using a three step procedure. Sub-sampling of the images were performed at factors 4,2,1 with knot-spacings 40, 20 and 10mm. Regularisation was set at λ = 0. 2, 0. 2 and 0. 4, and within each optimisation-step λ was a product of the tentative λ and the mean-sum-of-squared deviations, leading to a succesive relaxation of the regularisation. Running the method to converegence at full image resolution and a knot-spacing of 10mm reguires almost an hour on an intel Mac laptop. In this example the values of the Jacobian matrix was between 0. 3 and 3 for all voxels within all subjects, demonstrating that the inter-subjecr agreement was achieved without any "extreme" volume changes. We believe that this is largely attributed to the pyramid approach, and if attempting to register to the full resolution in a single step the warps tend to more "extreme".
section_name: RESULTS

section_content: We have described the implementation of a method/framework for small-displacement nonlinear registration of brain MR images. The results look promising and a first release of the Figure 2 : Transversal slice at z=56mm in the MNI-space of the same ten subjects as in figure 1. Average across affine registered images to the left and of non-linearly registered to the right. Note how the superior frontal sulcus, the precentral, central and postcentral sulci, the intraparietal sulcus and the cingulate gyrus are all clearly discernible in the non-linear average indicating a high degree of overlap. In the affine average in contrast it is really only the central and the cingulate gyrus that are discernible, and even those with large amount of blurring. This is further demonstrated by the manual traces of the central sulcus. 

software is scheduled to July 2007. Future work will aim at validation, finding the best set of parameters/schedules for different types of images and at a better modelling of the signal by including bias-field, and a physics based mapping of intensities between the two image (or between the template and the input image).
section_name: CONCLUSIONS

section_content: To see why this might all be useful let us look at a class derived from fnirt CF that implements the "sum of squared differences with scaling" cost-function described above.