The Objective Function 23

When calculating the correlation ratio, one should note that the width and the number of bins are important factors pertaining to the performance and the validity of CR. The perfor-mance of CR may drop significantly if the bin width is either too large or too small. Addi-tionally, the choice of the evaluation region is also important. To increase the efficiency, the evaluation region should be as compact as possible. Also, the calculation of CR(A, B, Φ) and CR(B, A, Φ⁻¹) involves different evaluation regions. For now we simply denote the evaluation region CR(A, B, Φ) as Ω and the evaluation region for CR(B, A, Φ⁻¹) as Ω^inv. Selection of the optimal bin width and the evaluation regions is further discussed in Chap-ter 4.

3.3.2 The Prior Term

The prior term is a measure which indicates how probable the transformation function is. This term can ensure the transformation to be realistic according to a certain prior knowledge such as smoothness of deformation or preservation of topology. Various types of prior terms have been used in previous works of registrations. This includes membrane energy, bending energy and linear-elastic energy [5].

The prior term used in our work is the membrane energy of the velocity field, which is also known as the Laplacian model:

E_Laplacian(V ) = 1 In here, Ω is the volume involved in the estimation. In our case, V is the support of the current RBF. The membrane energy is term that favors smoother deformation. This prior term was also used in other works [7, 22]. Unlike Liu et al. [22], in which the membrane energy was used to regularize the displacement field, this term was used to regularize the velocity field in the proposed algorithm. The user-specified weight (λ) controls how much the prior term effects the optimization. Large weight leads to results with smoother defor-mations and lower similarity measures. The registration results using small weight have

24 Methods

higher similarity but may be unrealistic. In our work, we set the value of λ empirically to 0.05 according to experiments using T1-MRIs.

3.4 Optimization

This section is divided into two parts. The first part explains how we designed the optimization problem as separate optimization problems in different local regions. The second part describes the symmetric optimization algorithm used in our work.

3.4.1 Local Optimization Scheme

In most registration approaches, coefficients of all basis functions are estimated simul-taneously. This involves optimization in an extremely high-dimensional parameter space.

Even with the use of basis functions to reduce the number of parameters, the number of pa-rameters to be estimated can still easily reach millions. The curse of dimensionality arises when finding the parameter set in such extremely high-dimensional search space, causing the typical high time complexity of registration algorithms. In the light of this problem, Rohde et al. [31] proposed a greedy approach which optimizes the basis functions sep-arately. In their work, 8 RBFs are applied to each identified region of mis-registration, and the vector coefficients of these 8 RBFs are estimated simultaneously. By doing so, the high-dimensional optimization problem is reduced to a sequence of 24-parameter op-timization problems. Liu et al. [22] further elevated this idea by optimizing all regularly deployed basis functions separately, which results in a sequence of 3-parameter optimiza-tion problems. The results of these works have showed that separating the optimizaoptimiza-tion of parameters yields significant improvement in term of the speed while maintaining decently high accuracy.

A similar greedy approach to Rohde et al. and Liu et al. is used in our work. In the proposed algorithm, a cumulative velocity field Vkis used to store the sum of all previously

3.4 Optimization 25

estimated RBFs from the start to the current step k:

V_k =

k

X

i=1

α_iρ_i(x) (3.19)

Each RBF is estimated such that the objective function will be minimized when this RBF is added to the cumulative velocity field:

α_i = arg min

αi

E(I_s, I_t, exp(V_i−1+ α_iρ_i(x))) (3.20) The estimated RBF is then added to the cumulative velocity field to form a new one:

V_i ← V_i−1+ α_iρ_i(x) (3.21)

Repeating these steps through all N RBFs yields the velocity field representing the overall transformation:

Φ = exp(V_N) (3.22)

Since RBFs are estimated sequentially, the order of RBFs in this optimization scheme may affect the accuracy of the final result. In our work, the order of RBFs is given according to the distance of each RBF from the brain center such that RBFs closer to the brain center are estimated before those farther from the brain center. This design is based on the fact that the cerebral cortex has higher structural complexity and anatomical variability than the structures near the center of the brain (e.g. corpus callosum). As a result, registering tissues near the brain surface is more difficult than registering tissues around the brain center. In the light of this fact, RBFs are designed to be estimated in the ascending order of their distances from the brain center. Registering around the brain center first may warp the cortex area to better initial positions, potentially leading to better estimation and thus higher accuracy.

3.4.2 Optimization Algorithm

The symmetry of the objective function alone still does not guarantee the symmetry of the whole registration algorithm. To make the registration algorithm symmetric, all elements in the algorithm must be unbiased toward the order within the image pair, i.e.,

26 Methods

which image in the image pair is the target image and which image is the source image.

Therefore, a symmetric optimization algorithm is necessary in our work.

The optimization algorithm used in our work is a modified version of the downhill sim-plex method. Downhill simsim-plex method [26] is an optimization algorithm that is also used by BIRT. This method is efficient but is rather unstable in terms of the initial positions of the simplex points. A slight change in the initial positions of simplex points may gener-ate a very different result. Furthermore, the original downhill simplex method initializes the simplex points in a random manner. As a result, BIRT cannot reproduce the same registration result in different trials using the same pair of image. This makes BIRT an unstable algorithm. To solve this problem, we use the gradient of the objective function at the origin of the orientation coefficient (α = [0 0 0]^T) to initialize the simplex points.

By doing so, we can ensure the symmetry of the optimization algorithm. An example is shown in Figure 3.7. According to the result of our expreiment, this symmetric downhill simplex method has similar speed and accuracy to the original downhill simplex method using random initial simplex points.