Synopsis

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.

INTRODUCTION

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"⁴ (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⁵. 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⁸ dataset.

METHOD DEVELOPMENT

The JSV penalty function (adapted from Ashburner⁷) 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⁹, 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¹⁰, and implemented as GPU code using CUDA¹¹. 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¹².

METHOD TESTING

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¹³) and nonlinear B-spline registration with Bending Energy¹⁴ (BE) regularisation (FNIRT) at 10mm and 1mm final knot spacings.

RESULTS

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.

DISCUSSION

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.

CONCLUSION

Acknowledgements

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.)

References

[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.