# Diffeomorphic Metric Mapping and Probabilistic Atlas Generation of Hybrid Diffusion Imaging based on BFOR Signal Basis

We propose a large deformation diffeomorphic metric mapping algorithm to align multiple b-value diffusion weighted imaging (mDWI) data, specifically acquired via hybrid diffusion imaging (HYDI), denoted as LDDMM-HYDI. We then propose a Bayesian model for estimating the white matter atlas from HYDIs. We adopt the work given in Hosseinbor et al. (2012) and represent the q-space diffusion signal with the Bessel Fourier orientation reconstruction (BFOR) signal basis. The BFOR framework provides the representation of mDWI in the q-space and thus reduces memory requirement. In addition, since the BFOR signal basis is orthonormal, the L2 norm that quantifies the differences in the q-space signals of any two mDWI datasets can be easily computed as the sum of the squared differences in the BFOR expansion coefficients. In this work, we show that the reorientation of the q-space signal due to spatial transformation can be easily defined on the BFOR signal basis. We incorporate the BFOR signal basis into the LDDMM framework and derive the gradient descent algorithm for LDDMM-HYDI with explicit orientation optimization. Additionally, we extend the previous Bayesian atlas estimation framework for scalar-valued images to HYDIs and derive the expectation-maximization algorithm for solving the HYDI atlas estimation problem. Using real HYDI datasets, we show the Bayesian model generates the white matter atlas with anatomical details. Moreover, we show that it is important to consider the variation of mDWI reorientation due to a small change in diffeomorphic transformation in the LDDMM-HYDI optimization and to incorporate the full information of HYDI for aligning mDWI.

## Authors

• 3 publications
• 3 publications
• 24 publications
• 2 publications
• 2 publications
• 3 publications
• 8 publications
10/10/2013

### Bayesian Estimation of White Matter Atlas from High Angular Resolution Diffusion Imaging

We present a Bayesian probabilistic model to estimate the brain white ma...
12/12/2018

### DeepTract: A Probabilistic Deep Learning Framework for White Matter Fiber Tractography

We present DeepTract, a deep-learning framework for estimation of white ...
07/24/2011

### Diffeomorphic Metric Mapping of High Angular Resolution Diffusion Imaging based on Riemannian Structure of Orientation Distribution Functions

In this paper, we propose a novel large deformation diffeomorphic regist...
06/03/2017

### Optimal Envelope Approximation in Fourier Basis with Applications in TV White Space

Lowpass envelope approximation of smooth continuous-variable signals are...
01/16/2016

### Estimation of Fiber Orientations Using Neighborhood Information

Data from diffusion magnetic resonance imaging (dMRI) can be used to rec...
05/19/2017

### Fiber Orientation Estimation Guided by a Deep Network

Diffusion magnetic resonance imaging (dMRI) is currently the only tool f...
09/25/2014

### Deconvolution of High-Dimensional Mixtures via Boosting, with Application to Diffusion-Weighted MRI of Human Brain

Diffusion-weighted magnetic resonance imaging (DWI) and fiber tractograp...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Diffusion-weighted MRI methods are promising tools for characterizing tissue microstructure. While diffusion tensor imaging (DTI) and high angular resolution diffusion imaging (HARDI) methods are widely used methods, they do not provide a complete description of the diffusion distribution. In order to more accurately reconstruct the ensemble average propagator (EAP), a thorough sampling of the

-space sampling is needed, which requires multiple -value diffusion weighted imaging (mDWI). The EAP estimation using mDWI better characterizes more complex neural fiber geometries and non-Gaussian diffusion behavior when compared to single -value techniques. Recently, new -space imaging techniques, diffusion spectrum imaging (DSI) Wedeen et al. (2005) and hybrid diffusion imaging (HYDI) Wu and Alexander (2007) have been developed for estimating the EAP. HYDI is a mDWI technique that samples the diffusion signal along concentric spherical shells in the

-space, with the number of encoding directions increased with each shell to increase the angular resolution with the level of diffusion weighting. Originally, HYDI employed the fast Fourier transform (FFT) to reconstruct the EAP. However, the recent advent of analytical EAP reconstruction schemes, which obtain closed-form expressions of the EAP, obviate the use of the FFT in HYDI. One such technique successfully validated on HYDI datasets is Bessel Fourier orientation reconstruction (BFOR)

Hosseinbor et al. (2013). While mDWI techniques like HYDI have not been widely used, the new human connectome project Essen et al. (2013) and spin-off projects will likely significantly increase the application. However, there is a lack of fundamental image analysis tools for mDWI and EAP maps, such as registration and atlas generation, that can fully utilize their information.

In the last decades, researchers have spent great efforts on developing registration algorithms to align diffusion tensors derived from DTI and orientation distribution functions (ODFs) derived from HARDI [e.g., Alexander et al. (2001); Raffelt et al. (2011); Du et al. (2012)]. However, registration algorithms directly based on DWIs are few. The direct alignment of DWIs in the -space utilizes the full diffusion information, is independent of the choice of diffusion models and their reconstruction algorithms (e.g., tensor, ODF), and unifies the transformation to align the local diffusion profiles defined at each voxel of two brains Dhollander et al. (2010); Yap and Shen (2012); Zhang et al. (2012). Dhollander et al.Dhollander et al. (2010) developed an algorithm that transforms the diffusion signals on a single shell of the -space and preserves anisotropic as well as isotropic volume fractions. Yap et.al Yap and Shen (2012) proposed to decompose the diffusion signals on a single shell of the -space into a series of weighted diffusion basis functions, reorient these functions independently based on a local affine transformation, and then recompose the reoriented functions to obtain the final transformed diffusion signals. This approach provides the representation of the diffusion signal and also explicitly models the isotropic component of the diffusion signals to avoid undesirable artifacts during the local affine transformation. Zhang et al. Zhang et al. (2012) developed a diffeomorphic registration algorithm for aligning DW signals on a single shell of the -space.

Only recently, Dhollander et al. Dhollander et al. (2011) aligned DWIs on multiple shells of the -space by first estimating transformation using a multi-channel diffeomorphic mapping algorithm, in which generalized fractional anisotrophy (GFA) images computed from each shell were used as mapping objects, and then applying the transformation to DWIs in each shell using the DWI reorientation method in Dhollander et al. (2010). This approach neglected possible influences of the DWI reorientation on the optimization of the spatial transformation. Hsu et al. Hsu et al. (2012) generalized the large deformation diffeomorphic metric image mapping algorithm Beg et al. (2005) to DWIs in multiple shells of the -space and considered the image domain and -space as the spatial domain where the diffeomorphic transformation is applied to. The authors claimed that the reorientation of DWIs is no longer needed as the transformation also incorporates the deformation due to the shape differences in the diffusion profiles in the -space. It is a robust registration approach with the explicit consideration of the large deformation in both the image domain and the -space. However, its computational complexity and memory requirement are high.

While limited research has been done for aligning the HYDI images, efforts on the white matter atlas from HYDI is even less. Only recently, Dhollander et al. Dhollander et al. (2011) employed their multi-channel diffeomorphic matching algorithm for the atlas generation using HYDI datasets. To our best knowledge, there is no research on probablistic atlas genertion for HYDI.

In this paper, we propose a new large deformation diffeomorphic metric mapping (LDDMM) algorithm to align HYDI datasets, denoted as LDDMM-HYDI, and then develop a Bayesian estimation framework for generating the brain atlas. In particular, we adopt the BFOR framework in representing the -space signal Hosseinbor et al. (2013). Unlike the diffeomorphic mapping of mDWIs in Hsu et al. Hsu et al. (2012), the BFOR signal basis provides the representation of the -space signal and thus reduces memory requirement. In addition, since the BFOR signal basis is orthonormal, the norm that quantifies the differences in the -space signals can be easily computed as the sum of the squared differences in the BFOR expansion coefficients. In this work, we will show that the reorientation of the -space signal due to spatial transformation can be easily defined on the BFOR signal basis. Unlike the work in Dhollander et al. (2011), we will incorporate the BFOR signal basis into the LDDMM framework and derive the gradient descent algorithm for solving the LDDMM-HYDI variational problem with explicit orientation optimization. Using this registration approach, we will further estimate the brain white matter atlas from the -space based on a Bayesian model. This probabilistic model is the extension of the previous Bayesian atlas estimation for scalar-based intensity images Ma et al. (2008). With the aids of the BFOR representation and reorientation of mDWIs introduced in this work, we show that it is feasible to adopt the previous Bayesian atlas estimation model for scalar-valued images Ma et al. (2008) to HYDI. As shown below, the main contributions of this paper are:

1. to seek large deformation for aligning HYDI datasets based on the BFOR representation of mDWI.

2. to derive the rotation-based reorientation of the -space signal via the BFOR signal basis. This is equivalent to applying Wigner matrix to the BFOR expansion coefficients, where Wigner matrix can be easily constructed by the rotation matrix (see Section 3.1).

3. to derive the gradient descent algorithm for the LDDMM-HYDI variational problem with the explicit orientation optimization. In particular, we provide a computationally efficient method for calculating the variation of Wigner matrix due to the small variation of the diffeomorphic transformation (see Section 3.4).

4. to show that the LDDMM-HYDI gradient descent algorithm does not involve the calculation of the BFOR signal bases and hence avoids the discretization in the -space.

5. to propose a Bayesian estimation model for the -space signals represented via the BFOR signal basis and derive an expectation-maximization algorithm for solving it (see Section 4).

## 2 Review: BFOR Signal Basis

According to the work in Hosseinbor et al. (2013), the -space diffusion signal, , can be represented as

 S(x,q) =Nb∑n=1NY∑j=1cnj(x)Ψnj(q) , (1)

where and respectively denote the image domain and -space. is the -th BFOR signal basis with its corresponding coefficient, , at . is given as

 Ψnj(q)=2√αnl(j)√πτ3Jl(j)+3/2(αnl(j))jl(j)(αnl(j)|q|τ)Yj(q|q|) . (2)

Here, is the root of the order spherical Bessel (SB) function of the first kind . is the radial distance in -space at which the Bessel function goes to zero. are the modified real and symmetric spherical harmonics (SH) bases as given in Goh et al. (2009). is the Bessel function of the first kind. is the number of terms in the modified SH bases of truncation order , while is the truncation order of radial basis. Note that in the original publication on BFOR Hosseinbor et al. (2013), the BFOR basis was not normalized to unity, but we have rectified this in Eq. (2). We refer readers to Hosseinbor et al. (2013) for more details on the BFOR algorithm.

Using the fact that the BFOR signal basis is orthonormal, the -norm of can be easily written as

 ∥S(x,q)∥2=√∫x∈R3∫q∈R3S2(x,q)dqdx= ⎷∫x∈R3Nb∑n=1NY∑j=1cnj(x)2dx . (3)

## 3 Large Deformation Diffeomorphic Metric Mapping for HYDI

### 3.1 Rotation-Based Reorientation of S(x,q)

We now discuss the reorientation of when rotation transformation is applied. We assume that the diffusion profile in each shell of the -space remains in the same shell after the reorientation. However, its angular profile in each shell of the -space is transformed according to the rotation transformation. Hence, we define

 RS(x,q)=S(x,|q|R−1q|q|) .

According to the BFOR representation of in Eq. (1), we thus have

 RS(x,q)=Nb∑n=1NY∑j=1cnj(x)2√αnl(j)√πτ3Jl(j)+3/2(αnl(j))jl(j)(αnl(j)|q|τ)Yj(R−1q|q|) .

This indicates that the rotation reorientation of mDWI is equivalent to applying the rotation transformation to the real spherical harmonics, . According to the work in Geng et al. (2011), the rotation of can be achieved by the rotation of their corresponding coefficients, yielding

 (4)

where is the th element of Wigner matrix constructed based on (see details in Geng et al. (2011)). We can see that the same Wigner matrix is applied to at a fixed . For the sake of simplicity, we rewrite Eq. (4) in the matrix form, i.e.,

 RS(x,q)=(M(R) c(x) )⊤Ψ(q) ,

where is a sparse matrix with diagonal blocks of .

is a vector that concatenates coefficients

in the order such that at a fixed , corresponds to . concatenates the BFOR signal basis.

### 3.2 Diffeomorphic Group Action on S(x,q)

We define an action of diffeomorphisms on , which takes into consideration of the reorientation in the -space as well as the transformation of the spatial volume in . Based on the rotation reorientation of in Eq. (4), for a given spatial location , the action of on can be defined as

 ϕ⋅S(x,q) =S(ϕ−1(x),R−1ϕ−1(x)q)

where can be defined in a way similar to the finite strain scheme used in DTI registration Alexander et al. (2001). That is, , where is the Jacobian matrix of at . For the remainder of this paper, we denote this as

 ϕ⋅S(x,q) =((M(Rx)c(x))⊤)∘ϕ−1(x)Ψ(q) , (5)

where indicates as the composition of diffeomorphisms.

### 3.3 Large Deformation Diffeomorphic Metric Mapping for HYDIs

The previous sections equip us with an appropriate representation of HYDI mDWI and its diffeomorphic action. Now, we state a variational problem for mapping HYDIs from one subject to another. We define this problem in the “large deformation” setting of Grenander’s group action approach for modeling shapes, that is, HYDI volumes are modeled by assuming that they can be generated from one to another via flows of diffeomorphisms

, which are solutions of ordinary differential equations

starting from the identity map . They are therefore characterized by time-dependent velocity vector fields . We define a metric distance between a HYDI volume of a subject and an atlas HYDI volume as the minimal length of curves in a shape space such that, at time , . Lengths of such curves are computed as the integrated norm of the vector field generating the transformation, where , where is a reproducing kernel Hilbert space with kernel and norm . To ensure solutions are diffeomorphic, must be a space of smooth vector fields. Using the duality isometry in Hilbert spaces, one can equivalently express the lengths in terms of , interpreted as momentum such that for each , , where we let denote the inner product between and , but also, with a slight abuse, the result of the natural pairing between and in cases where is singular (e.g., a measure). This identity is classically written as , where is referred to as the pullback operation on a vector measure, . Using the identity and the standard fact that energy-minimizing curves coincide with constant-speed length-minimizing curves, one can obtain the metric distance between the template and target volumes by minimizing such that at time . We associate this with the variational problem in the form of

 J(mt)= (6)

where is a positive scalar. quantifies the difference between the deformed atlas and the subject . Based on Eq. (3) and (5), is expressed in the form of

 E=∫x∈Ω∥∥(M(Rx) catlas(x))∘ϕ−1(x)−c(s)(x)∥∥22dx . (7)

### 3.4 Gradient of J with respect to mt

We now solve the optimization problem in Eq. (6) via a gradient descent method. The gradient of with respect to can be computed via studying a variation on such that the derivative of with respect to is expressed in function of . According to the general LDDMM framework derived in Du et al. (2011), we directly give the expression of the gradient of with respect to as

 ∇J(mt) = 2mt+ληt , (8)

where

 ηt=∇ϕ1E+∫1t[∂ϕs(kVms)]⊤(ηs+ms)ds , (9)

Eq. (9) can be solved backward given . is the partial derivative of with respect to .

In the following, we discuss the computation of . We consider a variation of as and denote the corresponding variation in as . Denote for the simplicity of notation. We have

 ∂E∂ϵ∣∣ϵ=0 =∫x∈Ω∂∥∥(M(Rϵx)catlas(x))∘(ϕϵ1)−1(x)−c(s)(x)∥∥22∂ϵ∣∣ϵ=0dx (10) =2∫x∈Ω⟨^c(x)∘ϕ−11−c(s)(x),∇⊤x^c(x)∘ϕ−11∂(ϕϵ1)−1∂ϵ∣∣ϵ=0⟩dxterm (A) +2∫x∈Ω⟨^c(x)∘ϕ−11−c(s)(x),(∂M(Rϵx)catlas(x)∂ϵ∣∣ϵ=0)∘ϕ−11⟩dxterm (B) .

As the calculation of Term (A) is straightforward, we directly give its expression, i.e.,

 Term (A)=−2∫x∈Ω⟨(Dxϕ1)−⊤∇x^c(x)(^c(x)−c(s)(ϕ1(x)))det(Dxϕ1),h⟩dx . (11)

This term is similar to that in the scalar image mapping case. It seeks the optimal spatial transformation in the gradient direction of image weighted by the difference between the atlas and subject’s images.

The computation of Term (B) involves the differential of with respect to rotation matrix and the variation of with respect to the small variation of . Let’s first compute the derivative of with respect to rotation matrix . According to the work in Cetingul et al. (2012), the analytical form of this derivative can be solved using the Euler angle representation of but is relatively complex. Here, we consider Wigner matrix and the coefficients of the BFOR signal basis together, which leads to a simple numeric approach for computing the derivative of with respect to rotation matrix , i.e., . Assume , where

is a skew-symmetric matrix parameterized by

. From this construction, is the tangent vector at on the manifold of rotation matrices and is also a rotation matrix. Based on Taylor expansion, we have the first order approximation of as

 M(~Rx)catlas(x)≈^c(x)+∇⊤Rx^c(x)δμ .

Hence, we can compute as follows. Assume to be skew-symmetric matrices respectively constructed from . We have

 (12)

It is worth noting that this formulation significantly reduces the computational cost for . Since is independent of spatial location , , , and are only calculated once and applied to all .

We now compute the variation of with respect to the small variation of . This has been referred as exact finite-strain differential that was solved in Dorst (2005) and applied to the DTI tensor-based registration in Yeo et al. (2009). Here, we directly adopt the result from Yeo et al. (2009) and obtain

 ∂Rϵx∂ϵ∣∣ϵ=0=−Fx3∑i=1[ri×(Dxh⊤)i] , (13)

where . denotes as the cross product of two vectors. denotes the th column of matrix . .

Given Eq. (12) and (13), we thus have

 Term (B) =−2∫x∈Ω⟨^c(x)∘ϕ−11−c(s)(x),(∇Rx^c⊤(x)Fx3∑i=1[ri×(Dxh⊤)i])∘ϕ−11⟩dx (14) =−2∫x∈Ωω⊤x3∑i=1[ri×(Dxh⊤)i]dx =−2∫x∈Ω3∑i=1⟨ωx×ri,∇xhi⟩dx ,

where

 ω⊤x=(∇Rx^c(^c(x)−c(s)(ϕ1(x))))⊤Fxdet(Dxϕ1) , (15)

and is approximated as

 Dxh=⎡⎢ ⎢⎣∇xh⊤1∇xh⊤2∇xh⊤3⎤⎥ ⎥⎦≈12Δd⎡⎢ ⎢⎣h1,xX+−h1,xX−h1,xY+−h1,xY−h1,xZ+−h1,xZ−h2,xX+−h1,xX−h2,xY+−h2,xY−h2,xZ+−h2,xZ−h3,xX+−h3,xX−h3,xY+−h3,xY−h3,xZ+−h3,xZ−⎤⎥ ⎥⎦,

where are the neighbors of in directions, respectively. is the distance of these neighbors to . Here, term (B) seeks the spatial transformation such that the local diffusion profiles of the atlas and subject’s HYDIs have to be aligned.

In summary, we have

 ∂E∂ϵ∣∣ϵ=0 ≈−2∫x∈Ω⟨(Dxϕ1)−⊤∇x^c(x)(^c(x)−c(s)(ϕ1(x)))det(Dxϕ1),h⟩dx (16) −1Δd∫x∈Ω3∑k=1⎧⎪ ⎪⎨⎪ ⎪⎩⟨ωx×rk,⎡⎢ ⎢⎣hk,xX+hk,xY+hk,xZ+⎤⎥ ⎥⎦⟩−⟨ωx×rk,⎡⎢ ⎢⎣hk,xX−hk,xY−hk,xZ−⎤⎥ ⎥⎦⟩⎫⎪ ⎪⎬⎪ ⎪⎭dx.

Therefore, can be obtained from Eq. (16).

### 3.5 Numerical Implementation

We so far derive and its gradient in the continuous setting. In this section, we elaborate the numerical implementation of our algorithm under the discrete setting. Since HYDI DW signals were represented using the orthonormal BFOR signal bases, both the computation of in Eq. (6) and the gradient computation in Eq. (16) do not explicitly involve the calculation . Hence, we do not need to discretize the -space. In the discretization of the image domain, we first represent the ambient space, , using a finite number of points on the image grid, . In this setting, we can assume to be the sum of Dirac measures, where is the momentum vector at and time . We use a conjugate gradient routine to perform the minimization of with respect to . We summarize steps required in each iteration during the minimization process below:

1. Use the forward Euler method to compute the trajectory based on the flow equation:

 dϕt(xi)dt = N∑j=1kV(ϕt(xi),ϕt(xj))αj(t) . (17)
2. Compute based on Eq. (16).

3. Solve in Eq. (9) using the backward Euler integration, where indices , with the initial condition .

5. Evaluate when , where is the adaptive step size determined by a golden section search.

## 4 Bayesian HYDI Atlas Estimation

The above registration problem requires an atlas. In this section, we introduce a general framework for Bayesian HYDI atlas estimation. Given observed HYDI datasets for , we assume that each of them can be estimated through an unknown atlas and a diffeomorphic transformation such that

 S(i)≈^S(i)=ϕ(i)⋅Satlas. (18)

The total variation of relative to is then denoted by . The goal here is to estimate the unknown atlas and the variation . To solve for the unknown atlas , we first introduce an ancillary “hyperatlas” , and assume that our atlas is generated from it via a diffeomorphic transformation of such that . We use the Bayesian strategy to estimate and from the set of observations by computing the maximum a posteriori (MAP) of . This can be achieved using the Expectation-Maximization algorithm by first computing the log-likelihood of the complete data () when are introduced as hidden variables. We denote this likelihood as . We consider that the paired information of individual observations, for , as independent and identically distributed. As a result, this log-likelihood can be written as

 logfσ(ϕ,ϕ(1),…,ϕ(n),S(1),…S(n)|S0) (19) = logf(ϕ|S0)+n∑i=1{logf(ϕ(i)|ϕ,S0)+logfσ(S(i)|ϕ(i),ϕ,S0)} ,

where

is the shape prior (probability distribution) of the atlas given the hyperatlas,

. is the distribution of random diffeomorphisms given the estimated atlas (). is the conditional likelihood of the HYDI data given its corresponding hidden variable and the estimated atlas (). In the remainder of this section, we first adopt and introduced in (Ma et al., 2008; Qiu et al., 2010) under the framework of LDDMM and then describe how to calculate in §4.2 based on the BFOR representation of HYDI.

### 4.1 The Shape Prior of the Atlas and the Distribution of Random Diffeomorphisms

In the framework of LDDMM-HYDI, we adopt the previous work in Ma et al. (2008); Qiu et al. (2010) and briefly describe the construction of the shape prior (probability distribution) of the atlas, . Under the LDDMM framework, we can compute the prior via , i.e.,

 f(ϕ|S0)=f(m0|S0) , (20)

where is initial momentum defined in the coordinates of such that it uniquely determines diffeomorphic geodesic flows from to the estimated atlas. When remains fixed, the space of the initial momentum provides a linear representation of the nonlinear diffeomorphic shape space, , in which linear statistical analysis can be applied. Hence, assuming is random, we immediately obtain a stochastic model for diffeomorphic transformations of . More precisely, we follow the work in (Ma et al., 2008; Qiu et al., 2010) and assume to be a centered Gaussian random field (GRF) model. The distribution of is characterized by its covariance bilinear form, defined by

 Γm0(v,w)=E[m0(v)m0(w)] ,

where are vector fields in the Hilbert space of with reproducing kernel .

We associate with . The “prior” of in this case is then

 1Zexp(−12⟨m0,kVm0⟩2),

where is the normalizing Gaussian constant. This leads to formally define the “log-prior” of to be

 logf(m0|S0)≈−12⟨m0,kVm0⟩2 , (21)

where we ignore the normalizing constant term .

We can construct the distribution of random diffeomorphisms, , in the similar manner. We define via the corresponding initial momentum defined in the coordinates of . We also assume that is random, and therefore, we again obtain a stochastic model for diffeomorphic transformations of . is assumed to be a centered GRF model with its covariance as , where is the reproducing kernel of the smooth vector field in a Hilbert space . Hence, we can define the log distribution of random diffeomorphisms as

 logf(ϕ(i)|ϕ,S0)≈−12⟨m(i)0,kπVm(i)0⟩2 . (22)

where as before, we ignore the normalizing constant term .

### 4.2 The Conditional Likelihood of the HYDI Data

Given the representation of the diffusion signals in the -space using BFOR bases, we construct the conditional likelihood of the HYDI data via the BFOR coeffcients. We assume that

has a multivariate Gaussian distribution with mean of

and covariance , where are the BFOR coefficients associated with .

From the Gaussian assumption, we can thus write the conditional “log-likelihood” of given and as

 logfσ(S(i)|ϕ(i)1,ϕ1,S0) (23) ≈

where as before, we ignore the normalizing Gaussian term.

### 4.3 Expectation-Maximization Algorithm

We have shown how to compute the log-likelihood shown in Eq. (19). In this section, we will show how we employ the Expectation-Maximization algorithm to estimate the atlas, , for , and . From the above discussion, we first rewrite the log-likelihood function of the complete data in Eq. (19) as

 logfσ(m0,m(1)0,…,m(n)0,S(1),…S(n)|S0) (24) ≈ −12⟨m0,kVm0⟩2−n∑i=1{12⟨m(i)0,kπVm(i)0⟩2

where and can be computed based on Eq. (5).

The E-Step. The E-step computes the expectation of the complete data log-likelihood given the previous atlas

and variance

. We denote this expectation as given in the equation below,

 Q(m0,σ2|mold0,σ2% old) (25) = E{logfσ(m0,m(1)0,…,m(n)0,S(1),…S(n)|S0)∣∣mold0,σ2old,S(1),⋯,S(n),S0} ≈ −12⟨m0,kVm0⟩2−n∑i=1E[12⟨m(i)0,kπVm(i)0⟩2 +

The M-Step. The M-step generates the new atlas by maximizing the -function with respcet to and . The update equation is given as

 mnew0,σ2new (26) = argmaxm0,σ2Q(m0,σ2|mold0,σ2old) = argminm0,σ2{⟨m0,kVm0⟩2

where we use the fact that the conditional expectation of is constant. We solve and by separating the procedure for updating using the current value of , and then optimizing using the updated value of .

We now derive how to update values of and from -function in Eq. (26). It is straightforward to obtain by taking the derivative of with respect to and setting it to zero. Hence, we have

 σ2new=1nn∑i=1∫x∈Ω∥∥∥(M(Rx)catlas(x))∘(ϕ(i)1)−1(x)−c(i)(x)∥∥∥2dx . (27)

For updating , let . Using the change of variables strategy, we have

 n∑i=1E[∫x∈Ω12σ2∥∥∥(M(Rx)catlas(x))∘(ϕ(i)1)−1(x)−c(i)(x)∥∥∥2dx] = n∑i=1E[∫y∈Ω12σ2∥∥∥M(Ry)catlas(y)−c(i)(y)∘ϕ1(y)∥∥∥2|Dϕ(i)1(y)|dy] = = +2⟨M(Ry)catlas(y)−¯¯¯c0(y),¯¯¯c0(y)−c(i)(y)∘ϕ1(y)⟩}|Dϕ(i)1(y)|]dy .

Since the second item in the above equation is independent of and ,we have

 mnew0=argminm0{⟨m0,kVm0⟩2+1σ2new∫y∈Ωα(y)∥∥∥(M