Frederik J. Lange^{1}, Stephen M. Smith^{1}, and Jesper L. R. Andersson^{1}

We present a method for regularising nonlinear registration which produces deformations which are biologically more plausible than conventional techniques. Our method, the Symmetric Prior for Regularisation of Elastic Deformations (SPRED), not only enforces diffeomorphism, but additionally penalises linear, planar and volumetric changes. Application of SPRED to the high quality NIREP dataset produced results whose quality matches that of established registration methods. The resulting deformations show significantly more plausible Jacobian distributions, both in terms of spatial locality and intensity. Future work will look to extend SPRED to include variable spatial priors, allowing different brain regions freedom to deform by varying amounts.

Many group-based analyses rely on nonlinear registration to transform an individual subjects' scans into a common space. When analysing the validity of the resulting transformation, an important aspect to consider is whether it is both "one-to-one" and "onto", meaning that all points map uniquely between spaces without folding or tearing^{1-3}. Such deformations are termed "diffeomorphic", and several "large-deformation frameworks"^{4} (theoretically) guarantee these conditions are met by generating deformations through the composition/integration of many small, smooth displacement fields^{2,3}.

Diffeomorphism alone should not, however, be considered sufficient to declare a transformation valid. Yet diffeomorphism is the only thing such techniques are capable of guaranteeing, as their formulation provides little or no direct control over, for example, the range of allowable volumetric changes. Highly biologically unlikely expansions or contractions are therefore possible, such as a functional group of cell bodies, or a bundle of white matter axons, being compressed by a factor ≥100.

Our work seeks to overcome this limitation by leveraging a "small-deformation framework"^{1,4}, wherein we explicitly penalise biologically implausible changes^{5}. Based on work by Ashburner et al^{6,7}, we use the Jacobian Singular Values (JSVs) of the deformation to penalise linear, planar and volumetric changes in a symmetric fashion (i.e. expansion and contraction are treated identically). Thus, not only do deformations remain diffeomorphic, but biologically plausible deformations are preferentially selected.

Here we present a novel, 3D, B-spline parametrised, Symmetric Prior for Regularisation of Elastic Deformations (SPRED). This has been developed as part of our new Multi-Modal Optimisation and Registration Framework (MMORF) tool. We demonstrate the performance of SPRED when applied to the NIREP^{8} dataset.

The JSV penalty function (adapted from Ashburner^{7}) implemented in SPRED is:

\begin{split}C\left(\vec{p}\right)&=\sum\limits_{xyz}c_i\left(\vec{p}\right)\\&=\sum\limits_{xyz}\lambda{v}\left(1+\left|\mathbf{J}_i\right|\right)\sum\limits_{n=1}^{3}\left(\log{s_n}\right)^2\\&\approx\sum\limits_{xyz}\lambda{v}\left(1+\left|\mathbf{J}_i\right|\right)\sum\limits_{n=1}^{3}\left.\left(s_n^2+\frac{1}{s_n^2}-2\right)\right/4\\&=\sum\limits_{xyz}\lambda{v}\left(1+\left|\mathbf{J}_i\right|\right)tr\left(\mathbf{J}_i^{T}\mathbf{J}_i+\left(\mathbf{J}_i^{-1}\right)^{T}\mathbf{J}_i^{-1}-2\mathbf{I}\right)/4\\\\\textrm{Where:}\\C&=\textrm{total
cost/penalty}\\\vec{p}&=\textrm{deformation
parameters}\\xyz&=\textrm{all voxels in
image}\\c_i&=\textrm{cost/penalty at voxel
}i\\\lambda&=\textrm{cost/penalty weighting}\\v&=\textrm{voxel
volume}\\\mathbf{J}_i&=\mathbf{J}_i\left(\vec{p}\right)=\textrm{Jacobian
matrix at voxel
}i\\s_n&=s_n\left(\vec{p}\right)=n^\textrm{th}\textrm{ singular
value of }\mathbf{J}_i\\\end{split}

This function conforms to our symmetric prior by equally penalising expansion and contraction of JSVs. Additionally it enforces diffeomorphism as the penalty tends to infinity when the JSVs approach either zero or infinity.

MMORF uses cubic B-splines to parametrise deformations and Gauss-Newton (GN) for optimisation. As this function is not of the $$$y=g\left(x\right)^2$$$ form required for GN assumptions to hold^{9}, a variable substitution is needed which, together with the chain rule, leads to the following gradient and Hessian formulae:

\begin{split}C\left(\vec{p}\right)&=\sum\limits_{xyz}g_i^2\left(\vec{p}\right)=\sum\limits_{xyz}c_i\left(\vec{p}\right)\\\\\nabla_C\left(\vec{p}\right)&=2\sum\limits_{xyz}\frac{\partial{g_i}}{\partial{\vec{p}}}g_i\left(\vec{p}\right)=\sum\limits_{xyz}\frac{\partial{c_i}}{\partial{\vec{p}}}\\&=\frac{\partial{c_i}}{\partial{\vec{s}}}\times\frac{\partial{\vec{s}}}{\partial{\vec{j}}}\times\frac{\partial{\vec{j}}}{\partial{\vec{p}}}\\\\H_C\left(\vec{p}\right)&=2\sum\limits_{xyz}\frac{\partial{g_i}}{\partial{\vec{p}}}\left(\frac{\partial{g_i}}{\partial{\vec{p}}}\right)^T=\sum\limits_{xyz}\frac{1}{2c_i\left(\vec{p}\right)}\frac{\partial{c_i}}{\partial{\vec{p}}}\left(\frac{\partial{c_i}}{\partial{\vec{p}}}\right)^T\\&=\sum\limits_{xyz}\frac{1}{2c_i\left(\vec{p}\right)}\left(\frac{\partial{c_i}}{\partial{\vec{s}}}\times\frac{\partial{\vec{s}}}{\partial{\vec{j}}}\times\frac{\partial{\vec{j}}}{\partial{\vec{p}}}\right)\left(\frac{\partial{c_i}}{\partial{\vec{s}}}\times\frac{\partial{\vec{s}}}{\partial{\vec{j}}}\times\frac{\partial{\vec{j}}}{\partial{\vec{p}}}\right)^T\\\\\textrm{Where:}\\\vec{s}&=\textrm{vectorised
singular values of }\mathbf{J}_i\\\vec{j}&=\textrm{vectorised
elements of }\mathbf{J}_i\\\end{split}

The term $$$\frac{\partial{c_i}}{\partial{\vec{s}}}\times\frac{\partial{\vec{s}}}{\partial{\vec{j}}}$$$ is independent for each voxel making this calculation trivially parallelisable. A closed form solution for this term was derived using the Matlab Symbolic Toolbox^{10}, and implemented as GPU code using CUDA^{11}. The term $$$\frac{\partial{\vec{j}}}{\partial{\vec{p}}}$$$ is also calculated on the GPU using the method underpinning our previous work on movement by susceptibility distortion correction^{12}.

MMORF with SPRED was used to register all 240 possible combinations of the 16 subject NIREP dataset. A 6-stage, multi-resolution, multi-smoothing approach was used to a final knot spacing of 5mm.

The resulting transformations were then applied to 32 segmented grey matter regions for each combination. Jaccard coefficients of the overlap between transformed and reference regions were calculated.

The procedure was repeated using linear registration (FLIRT^{13}) and nonlinear B-spline registration with Bending Energy^{14} (BE) regularisation (FNIRT) at 10mm and 1mm final knot spacings.

All 240 SPRED regularised transformations remained entirely diffeomorphic. Despite attempting to project the final deformation onto the closest diffeomorphic representation, FNIRT still occasionally produced voxels with negative Jacobian determinants.

The Jaccard coefficients (combined left and right hemisphere) in Figure 1 show a consistent trend across regions, with nonlinear registration outperforming linear as expected. Figure 2 visually confirms that the registration was successful.

Figure 3 and 4 show that despite the data driving the registration in a similar direction, there are distinct differences is the Jacobian determinant distributions between SPRED and BE regularisation.

The Jaccard coefficients demonstrate that MMORF's performance and variability are comparable to the well established FNIRT, indicating that the strict SPRED constraint has not detrimentally affected performance. High visual similarity and low residual difference between registered images confirms this.

Figures 1 and 2 clearly show that SPRED produces high quality registration and figures 3 and 4 show that it does so with a much tighter distribution of Jacobian determinants than when using BE. We therefore posit that the SPRED registrations are more biologically plausible.

The Wellcome Centre for Integrative Neuroimaging is supported by core funding from the Wellcome Trust (203139/Z/16/Z).

We also gratefully acknowledge the support from the Wellcome-Trust Strategic Award 098369/Z/12/Z (J.L.R.A.) and the NIH Human Connectome Project (1U01MH109589-01 and 1U01AG052564-01, J.L.R.A.)

Additionally this work was supported by funding from the Oppenheimer Memorial Trust, St Catherine's College Oxford, The Rotary Foundation, and the Engineering and Physical Sciences Research Council (EPSRC) and Medical Research Council (MRC) [grant number EP/L016052/1] (F.J.L.)

[1] J. L. R. Andersson, M. Jenkinson, and S. Smith, “Non-linear registration, aka spatial normalisation. FMRIB Technial Report TR07JA2.,” 2007.

[2] M. F. Beg, M. I. Miller, A. Trouvé, and L. Younes, “Computing Large Deformation Metric Mappings via Geodesic Flows of Diffeomorphisms,” Int. J. Comput. Vis., vol. 61, no. 2, pp. 139–157, Feb. 2005.

[3] J. Ashburner, “A fast diffeomorphic image registration algorithm,” Neuroimage, vol. 38, no. 1, pp. 95–113, 2007.

[4] M. Miller et al., “Statistical methods in computational anatomy.,” Stat. Methods Med. Res., vol. 6, no. 3, pp. 267–299, 1997.

[5] T. Rohlfing, C. R. Maurer, D. A. Bluemke, and M. A. Jacobs, “Volume-Preserving Nonrigid Registration of MR Breast Images Using Free-Form Deformation With an Incompressibility Constraint,” vol. 22, no. 6, pp. 730–741, 2003.

[6] J. Ashburner, J. L. Andersson, and K. J. Friston, “High-dimensional image registration using symmetric priors.,” Neuroimage, vol. 9, no. 6 Pt 1, pp. 619–28, 1999.

[7] J. Ashburner, J. L. Andersson, and K. J. Friston, “Image registration using a symmetric prior--in three dimensions,” Hum. Brain Mapp., vol. 9, no. 4, pp. 212–225, 2000.

[8] G. E. Christensen et al., “Introduction to the Non-rigid Image Registration Evaluation Project (NIREP),” Lect. Notes Comput. Sci. (including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics), vol. 4057 LNCS, pp. 128–135, 2006.

[9] P. Chen, “Hessian Matrix vs. Gauss-Newton Hessian Matrix,” Siam, vol. 49, no. 4, pp. 1417–1435, 2011.

[10] Mathworks, “Symbolic Math Toolbox User’s Guide V8.1.” 2018.

[11] NVIDIA, “CUDA C Programming Guide.” 2017.

[12] F. Lange, M. Graham, I. Drobnjak, H. Zhang, J. Campbell, J. Andersson, "Estimating susceptibility-induced field changes directly from diffusion MRI imagesand overcoming associated computational bottlenecks through GPU parallelisation", ISMRM 2018.

[13] M. Jenkinson, C. F. Beckmann, T. E. J. Behrens, M. W. Woolrich, and S. M. Smith, “FSL,” Neuroimage, vol. 62, no. 2, pp. 782–790, 2012.

[14] F. Bookstein, “Quadratic variation of deformations,” Inf. Process. Med. Imaging, pp. 15–28, 1997.

Jaccard coefficients for the 16 combined left and right grey matter segmented regions of the NIREP dataset. All 240 combinations of registrations for the 16 subjects were used. We compare the results between one linear (FLIRT) and three nonlinear (FNIRT at 10mm and 1mm warp resolution, and MMORF at 5mm warp resolution) implementations. FNIRT employed a conventional Bending Energy regularisation whilst MMORF used our novel SPRED penalty. MMORF compares favourably against FNIRT both in terms of average performance and variance, but has the benefit of intrinsically limiting the Jacobian determinants to be diffeomorphic.

An example of a moving image (b) transformed to match a reference image (d). Visually, linear registration (c) does not greatly improve the perceived similarity, whereas after nonlinear registration using MMORF with SPRED (e) the similarity is markedly improved. The change in the residual error before (a) and after (f) nonlinear registration confirms the observed subjective improvement in alignment. The nonlinear result does not appear unrealistically over-deformed which suggests SPRED is sufficiently regularising the deformation.

Log-Jacobian determinant maps for conventional Bending Energy regularisation (a) and for SPRED regularisation (b). Both registrations were carried out to 5mm warp resolution using our MMORF registration tool until the Sum-of-Squared Differences costs were approximately equal. Bending Energy regularisation resulted in approximately 10,000 voxels having negative Jacobian determinants. This is in contrast to SPRED where all Jacobian determinants remained positive. Visually it is clear that, over and above remaining diffeomorphic, both the numerical range and the spatial extent of more extreme Jacobian determinants is reduced under SPRED regularisation.

Probability density for the values of the Jacobian determinant in the example registration from Figure 2. Here it can be seen that some of the Jacobian determinants are negative under Bending Energy regularisation, whereas no Jacobian determinants are negative under SPRED regularisation. Note that SPRED also displays a much narrower distribution of Jacobian determinants, indicating that deformations which require less volumetric changes are being preferentially selected during the registration process.