A Statistical Model for Simultaneous Template Estimation, Bias Correction, and Registration of 3D Brain Images

TL;DR

This paper introduces a generative statistical model for brain image analysis that simultaneously estimates templates, corrects bias fields, and registers images, improving robustness to variability in MRI data.

Contribution

It presents a novel joint inference framework for template estimation, bias correction, and registration, addressing limitations of existing methods that rely on bias-invariant measures.

Findings

Model accurately captures intensity heterogeneity in MRI images.

Provides reliable template estimation from registration.

Demonstrates effectiveness on synthetic and real data.

Abstract

Template estimation plays a crucial role in computational anatomy since it provides reference frames for performing statistical analysis of the underlying anatomical population variability. While building models for template estimation, variability in sites and image acquisition protocols need to be accounted for. To account for such variability, we propose a generative template estimation model that makes simultaneous inference of both bias fields in individual images, deformations for image registration, and variance hyperparameters. In contrast, existing maximum a posterori based methods need to rely on either bias-invariant similarity measures or robust image normalization. Results on synthetic and real brain MRI images demonstrate the capability of the model to capture heterogeneity in intensities and provide a reliable template estimation from registration.

Equations

$$I_i=\theta (\text{Exp}(v(w_i))+x_i+{ϵ}_i$$

$$I_i\approx {\theta }^{w_i^0}+Z_i(w_i-w_i^0)+x_i+{ϵ}_i,\text{where,}$$

$${\theta }^{w_i^0}=\theta (\text{Exp}(v(w_i^0))),Z_i={\nabla }_x{\theta (\text{Exp}(v(w_i^0)))}^T{|}_{x=\text{Exp}(v(w_i^0))}{J}_w\text{Exp}(v(w_i^0))\in {\mathbb{R}}^{n\times n_w},$$

$$w_i\sim \mathcal{N}(0,{\sigma }^2{C}_i),C_i={\mathbb{I}}_i\otimes C_0,x_i\sim \mathcal{N}(0,{\sigma }^2{S}_i),S_i={\mathbb{I}}_i\otimes S_0,$$

$${ϵ}_i\sim \mathcal{N}(0,{\sigma }^2{\mathbb{I}}_i).$$

$$\begin{array}{rl}& \mathcal{L}({\sigma }^2,S,C)=nk\text{log}{\sigma }^2+\sum _{i=1}^k\text{log}\text{det}{V}_i+\\ & \frac{1}{{\sigma }^2}\sum _{i=1}^k(I_i-{\theta }^{w_0}+Z_i{w}^0{)}^T{V}_i^{-1}(I_i-{\theta }^{w_0}+Z_i{w}^0).\end{array}$$

$$\begin{array}{rl}& \mathcal{L}({\sigma }^2,S,C)\approx nk\text{log}{\sigma }^2+\sum _{i=1}^k\sum _{j=1}^p\text{log}\text{det}{V}_{i,j}+\\ & \frac{1}{{\sigma }^2}\sum _{i=1}^k\sum _{j=1}^p(I_{i,j}-{\theta }^{w_0}+Z_{i,j}{w}^0{)}^T{V}_{i,j}^{-1}(I_{i,j}-{\theta }_j^{w_0}+Z_{i,j}{w}^0).\end{array}$$

$$\mathcal{P}(w)=\sum _{i=1}^k\sum _{j=1}^p(I_{i,j}-{\theta }_j^{w_i}{)}^T(S_j+I_j{)}^{-1}(I_{i,j}-{\theta }_j^{w_i})+w_i^T{C}_i{w}_i.$$

$$E[x_i|w_i^0,I_i]=S(S+I{)}^{-1}(I_i-{\theta }_j^{w_i}).$$

$$\begin{array}{rl}\frac{\partial \varphi (x,t)}{\partial t}& =v(\varphi (x,t)),\phi =\varphi (x,1)=\text{Exp}(v).\end{array}$$

$$v(x)=\sum _{m=1}^R{v}_m=\sum _{m=1}^R\sum _i^{N_m}{K}_m(c_i^m,x)w_i^m.$$

Full text

1 Introduction

Brain template estimation is becoming increasingly important since it facilitates a variety of applications such as segmentation, registration, or providing common coordinate systems for statistical analysis of shape models for a given population. At the core of template estimation methods is image registration. In a statistical setting, conventional image registration methods are often approached from a Bayesian viewpoint where one maximizes the posterior given the image data and a regularizing prior [1]. To avoid choosing parameters for controlling the deformation (regularization) in an ad-hoc fashion, recent methods have employed Bayesian models where parameters are estimated in a data-driven fashion by treating them as latent random variables drawn from a distribution with a smooth covariance structure [2, 3]. These methods employ $L^{2}$ similarity measures that are fragile towards deviations in model assumption; for instance, bias fields. To achieve reasonable results under model deviations, strong penalization on the variation in deformation parameters can be imposed, or similarity measures that are invariant to bias fields (e.g. mutual information) can be used instead of the $L^{2}$ data terms.

In this paper, we propose a statistical mixed-effects model where deformation and spatially correlated variation in intensity are modeled as random effects. This way, effects from deformations and variations in intensities due to scanners can be separately handled, and hyperparameters can be estimated in a data-driven fashion using maximum likelihood. We perform simultaneous estimation and prediction in the model to avoid bias in the estimation that can result from treating warping as a preprocessing step [4]. In addition, we propose a different estimation procedure in comparison to existing Bayesian methods. While current methods use non-linear sampling for marginalizing over the latent variables, we propose to use successive linearization around predictions of the latent variables allowing estimation with linear mixed-effect theory. This paper is built on the methods proposed by Raket et. al, [4] for analysis of 1D functional data.

Current probabilistic template estimation methods do not model multi-scale behavior of the deformations. Although innumerable multi-scale deformation models have been proposed for image registration, their application in template estimation still needs maturity. To this end, we propose to model deformations at different scales as latent variables that are drawn from different distributions with different covariance structure. Concretely, we utilize the kernel bundle structure [5, 6] to model velocity fields, and the multi-scale nature is interleaved in the covariance structure of the distribution of parameters at each kernel bundle level.

The contributions of the paper are as follows:

• •

We demonstrate the utility of a computationally feasible class of non-linear mixed effect models in 3D brain template estimation.

• •

We propose a model that handles effects from deformation and intensity variations separately and allows simultaneous estimation of template and variance hyperparameters and prediction of deformation and bias fields.

• •

We propose an iterative linearization of the model in the non-linear random effects that enables efficient maximum likelihood estimation of variance parameters.

• •

Finally, we propose the incorporation of scales in the deformation distribution via the kernel bundle representation.

2 Background

A number of registration-based approaches to template estimation have been proposed where the central aspect addressed was the choice of target co-ordinate system. Initial approaches, which include the popular minimum deformation template [7], choose a random image as a target co-ordinate system, and is iteratively updated by registering other images to it. This approach has been shown to be significantly biased towards the choice of the random image [8].

As an alternate approach, several papers have proposed the strategy of registering several images to a template, which is simultaneously estimated via an alternating optimization scheme [9, 10]. In a probabilistic formulation, the image matching term can be viewed as a log likelihood and the regularization as a log prior. In this direction, recent methods like [3] have employed tools such as non-linear mixed effects models to deal with the population effects (template or the fixed effect) and individual effects (deformations or the random effect). Methods have moved towards simultaneously modeling the template and inferring parameters of the deformation such as the regularization factor with expectation-maximization [2]. All the aforementioned methods rely on bias-invariant similarity measures or robust image normalization methods. To address intensity bias factors, Hromatka et. al, [11]., proposed to use a hierarchical Bayesian model for template estimation where two transformations are concatenated; one taking an individual image to an atlas of a site and the other that takes this warped image to the global atlas. However, this model takes into consideration very little about intra-site bias variations.

3 Statistical model

Consider a population of images $I_{i}:\mathbb{R}^{3}\rightarrow\mathbb{R},i=1\dots k$ and let $\theta$ be a template of these images; both measured on a discrete grid $\Omega\in\mathbb{Z}^{3}$. The individual observed image may be then defined in terms of the fixed and random effects as

[TABLE]

where the template $\theta$ is the fixed effect. More control over the smoothness of $\theta$ may be incorporated by specifying a parametric subspace for $\theta$. However, such constructions will not be entertained in this paper. The remaining effects are all random: The deformation $\text{Exp}(v(w_{i}))$ (Section 3.2) is a random field controlled by latent random deformation parameters $w_{i}$. A key contribution of this paper is the incorporation of the random spatially correlated effect $x_{i}$ that models a bias field. Note that $x_{i}$ is defined in the frame of the individual image and not the template. Converse constructions are also possible. Finally, $\epsilon_{i}$ is the i.i.d noise.

Following [4] and assuming that $\theta$ is smooth, we can linearize (1) around deformation parameters $w_{i}^{0}$ resulting in the linear model

[TABLE]

Here $\sigma^{2}$ is the noise variance, and the spatial covariance of the deformations and bias fields are controlled by the matrices $C$ and $S$. Note that $\nabla_{\mathbb{x}}$ denotes derivatives with respect to spatial coordinates. Also, note that deformations in the linear model are parameters of $w$ while the linearization point $w^{0}$ is iteratively optimized for during the estimation process. The model here is a single scale model where the covariances are constructed with only one scale of the kernel. Multi-scale version of this will be discussed in Section 3.2.

The first step of the analysis is to estimate $\theta$ with an initial guess of linearization point $w$. $\theta$ is computed by back-warping the images with the velocity field as $\theta^{w}_{0}=\frac{1}{k}\sum_{i=1}^{k}I_{i}(\text{Exp}(-v(w_{0}))).$ This simplified formulation of the conditional maximum-likelihood is a result of Henderson’s mixed-model [12] that simplifies when all observations are on a common grid. With this estimate of $\theta$, we estimate the variance parameters by minimizing the double negative log-likelihood of the linearized model

[TABLE]

where $V_{i}=Z_{i}^{T}C_{i}Z_{i}+S+\mathbb{I}_{i}$. Computing (2) directly is computationally intractable. This is because the dimensionality of $S,~{}V$ is $m^{2}$, where $m$ is the number of voxels in the image. We handle this by two assumptions: a) The support of the kernel used to construct the covariance matrix of the spatial correlation effect is sufficiently large that the inverse of the covariance matrix resembles a block-diagonal matrix, see e.g. assumption [13]; and b) the interaction between the blocks is negligible. This allows the likelihood to be approximated by an integral over smaller computationally tractable patches over the image. Therefore, (2) can be rewritten as:

[TABLE]

where $V_{i,j}=Z_{i,j}^{T}C_{i}Z_{i,j}+S_{j}+\mathbb{I}_{j}$. Here $j=1\dots p$ is the patch number of the image with total $p$ patches. The estimation process is now repeated: Given the current estimate of the template $\theta$ and the variance estimates, the new deformation parameters $w^{0}$ for the linearization point are chosen as the most likely predictions in the original non-linear model (1). That is, they are given by minimizing the negative log posterior of the deformation given the image data in model (1):

[TABLE]

To reconstruct the bias field for a given iteration, the conditional expectation of the spatially correlated effect $x_{i}$ given the data $I_{i}$ and the most likely deformation parameters $w_{i}^{0}$ are computed under the maximum likelihood estimates for the model parameters:

[TABLE]

Note that the best linear unbiased predictor (BLUP) [14] of $w_{i}$ in the linearized model given the image data $I_{i}$ is realized by $E[w_{i}|I_{i}]$ which is computed as $(C_{i}^{-1}+Z_{i}^{T}(\mathbb{I}_{i}+S)^{-1}Z_{i})^{-1}Z^{T}(\mathbb{I}_{i}+S)^{-1}(I_{i}-\theta^{w_{i}^{0}}+Zw^{i}_{0})$.

3.1 Covariance matrices

A key aspect of this statistical model is the choice of covariance matrices for the spatially correlated effects and the deformation effects. Current Bayesian methods model the inverse of the covariance matrix directly by the means of an operator; typically of the Cauchy-Navier type which takes the form, $L=-\alpha\nabla+\beta$ where $\nabla$ is the Laplace operator. The parameter $\alpha$ here controls the smoothness of the covariance. In this paper, since we work with patches, it is more intuitive to model the covariance matrix directly. We model the covariance matrices for both the deformation parameters and the spatial correlation using Wendland kernels [15] that are compactly supported reproducing kernels. The covariance matrices are constructed using $C=\lambda^{2}K(c_{i},c_{j})\ ,\ \lambda\in(0,\infty)$, $K(a,b)=r^{4}(4r+1),~{}r=(1-t,0)_{+},~{}t=\frac{||a-b||}{s}$, where $c_{i}$ are the kernel centers.

For the spatially correlated effect ($S=\beta K$), a similar representation is chosen. However, we choose the parametric subspace to be the same size as that of the image i.e., a kernel is centered at every voxel as opposed to the much smaller subspace spanned by the deformation kernels. The amplitude of these covariance matrices are controlled by parameters $\lambda,\beta$. These parameters are estimated by optimizing the likelihood (2). The smoothness of the deformation is controlled by $C$. We here keep $C$ fixed though parameters of $C$, e.g. range and scale, can also be optimized for in (2).

3.2 Multi-scale deformation model

We model deformation fields as path of diffeomorphisms generated by integrating stationary velocity fields (SVFs). Let $G\subset\text{Diff}(\Omega)$ be a Lie subgroup of the group of diffeomorphic transformations $\varphi:\Omega\rightarrow\Omega$, and let $\phi:\Omega\times\mathbb{R}\rightarrow\Omega$ be a path in $G$. Let $V$ be the tangent space of $G$ at identity Id containing velocity fields ${v}:\Omega\rightarrow\mathbb{R}^{d}$. In SVFs, the governing differential equation can be written as:

[TABLE]

The final transformation $\phi$ is the Lie group exponential $\text{Exp}({v})$. The integration in (6) can be numerically realized as an Euler integration. We use the kernel bundle framework [5, 6] in this paper. In short, the concept of the space of velocity fields $V$ is extended to a family $\hat{V}$ of spaces of velocity fields. The velocity fields are linear sums of individual kernels at $R$ levels

[TABLE]

Here $K$ is an interpolating Wendland kernel, $R$ is the total number of kernel bundle levels, $N_{m}$ is the number of kernels at each level and $c$ is center of each kernel at the kernel bundle level. The parameter $w^{m}$ is associated with the $m$th kernel bundle level, and we assume that $w_{i}^{m}\sim\mathcal{N}(0,\sigma^{2}C^{m})$, where $C^{m}$ is the covariance matrix for each kernel bundle level with its distinct support and smoothness.

4 Experiments and Results

We perform an evaluation on the MGH10 dataset111www.mindboggle.info. The dataset contains 10 images each with the dimension of 182x218x182 and a voxel resolution of 1,1.33,1 mm. The images are initially co-registered rigidly. For bias recovery only a subset of 5 subjects is used.

One of the major challenges in the field of image template estimation is validation. Since the underlying geometry of the space of images is unknown, the definition of modes of populations is difficult. Therefore, the validation needs to be tied to specific applications of template estimation like image segmentation or even validation of the underlying image registrations. For the 2 experiments in this paper, we chose 3 kernel bundle levels each with a control point spacing of 20, 10, 6 mm respectively. The support of the kernel used to construct the covariance matrix for the deformations is fixed to 4 across the kernel bundle levels. The support of the spatial correlation variance kernel was set to 40. We perform 2 experiments to validate our template estimation method.

4.1 Bias field recovery

We select a subset of 5 images from the database and add artificial multiplicative bias to the image of form $\text{bias}=1+\text{exp}(-\frac{\mathbb{x}^{2}+\mathbb{y}^{2}+\mathbb{z}^{2}}{2*30^{2}})*0.05.$. We then apply our statistical model to estimate the template and inspect whether the corrupted image can be recovered. As illustrated in Figure 1 the image recovered is free of the bias field illustrating the robustness of the method towards scanner-related artifacts.

4.2 Template estimation

In this paper, we evaluate the overlaps of segmentations of images mapped to the estimated template. This experiment was chosen to illustrate the performance of the underlying registration method in the paper. Figure 2 illustrates the mean and target overlaps estimated by mapping the 10 individual images to the template space and measuring pairwise overlaps. The overlaps are comparable to what state-of-the-art registration methods have obtained on the same dataset [16]. The kernel bundle scales are sequentially optimized. However, if one switches to a parallel optimization across scales, a significant improvement in the overlaps may be expected. Also in the illustration is the non-aligned template image and aligned template image. The sharpness of the template is particularly visible in corpus collosum.

To further demonstrate the effectiveness of the proposed template estimation method and benefits of multi-scale deformation, we visually inspect the atlas estimated in key regions like hippocampus and putamen. As illustrated in Figure 2, as the kernel bundle resolution becomes finer, the sharper are the boundaries of the anatomical structure.

5 Discussion and Conclusion

We presented a class of non-linear mixed effect models for estimating 3D brain templates, and we proposed to simultaneously estimate both bias field and deformation parameters in a data-driven maximum likelihood setting. In addition, we proposed to incorporate the kernel bundle framework in the random deformation effects to account for deformations occurring at different scales. We illustrated the application of our in both bias field correction and template estimation. The basis for the simultaneous estimation of the deformation and spatially correlated bias was results of [4] where the authors showed that estimation of deformation parameters is biased if the data is preprocessed first followed by the prediction of deformation parameters. We modeled deformations as random geodesics on the diffeomorphism group and estimated variance parameters of both the deformation and bias field from data instead of setting them ad-hoc. In future work, we will investigate replacing random geodesics by more geometrically natural distributions on the space of diffeomorphism [17]. An important next step for this work will be to extend the method to account for multiple population means.