Optimization

The voxel similarity measures described in this chapter all need to be incorporated into an iterative optimization scheme in order to register images. Depending on the similarity measure being used, it may or may not be possible to calculate derivatives of the similarity measure with respect to the parameters of the spatial transformation model (degree of freedom). If such derivatives can be calculated, then a calculus-based optimization algorithm such as those described elsewhere in this book can be used. Alternatively, non-calculus-based techniques are necessary. This section concentrates on optimizing information theoretic similarity measures that have been the focus of this chapter.

5.1 Discrete and Continuous Formulations of Mutual Information

There are two distinct formulations of mutual information registration. These methods differ in the way that they estimate the joint probability distribution function. The method described earlier calculates the joint PDF from the discrete joint histogram. The alternative method proposed by Viola and Wells [1,20], involves using a Parzen window density estimate, which generates a continuous PDF. Two randomly chosen samples of the voxels in the image pair are chosen. The first sample a of Na voxels is used to estimate the density function, and the second sample ft of Ng is used to estimate the entropy. The Parzen window function chosen by Wells ef aZ.  is G^, a

Gaussian with variance The estimated entropy of a distribution z is therefore

This approach can be used to calculate both the joint entropy and the marginal entropies. An advantage of this approach compared to the histogram method of PDF estimation is that it is possible to calculate the derivatives of the entropies with respect to the transformation parameters of T. A further advantage of this approach is that the PDF can be calculated using a small number of samples from the images. Wells suggests Na = Nf = 50.

5.2 Image Resampling

In order to achieve high registration accuracy by maximization of mutual information, it is important to calculate changes in the value of mutual information for small changes in the transformation parameters. It is therefore essential that the image resampling required for each transformation estimate not lead to artifacts in the PDF that would result in local extrema in the search space. Maes et al.  propose the use of partial volume interpolation as a means of accurately estimated the joint histogram, and hence PDF. In general, an image transformation will require interpolation between existing sample points, and fast interpolation algorithms such as nearest-neighbor and trilinear interpolation are well known to introduce errors . In partial volume interpolation, the contribution to the joint histogram from a voxel in image A and the corresponding location in the transformed image B is not a single value calculated by trilinear interpolation of image B, but multiple values distributed over all eight nearest-neighbor intensity values, using the same weights as used in trilinear interpolation. This approach can introduce other artifacts and has not been universally adopted. Many authors [21,23-25] find trilinear interpolation sufficient.

5.3 Capture Range

In many applications of optimization algorithms, the goal is to find the global optimum value in an infinite parameter space, and a successful algorithm will manage to avoid finding solutions in local optima or on plateaus in the search space. This tends to be the situation for image registration using surface matching, but for image registration using voxel similarity measures, the situation is quite different. The similarity measures described in this chapter give an indication of how well the data is matched within the volume of overlap of the images. If the amount of overlap is small, then the measures can give a misleading indication of registration accuracy. For example, taking typical images in which there is air at the edge of the field of view in both modalities, if these are misaligned so much that the only overlapping voxels contain air in both modalities, then voxel similarity measures will judge this to be a good alignment. In fact, it is common for such an alignment to be globally optimal as far as the similarity measure is concerned, but clearly incorrect. A successful voxel similarity measure will have a limited range of transformations for which the value of the measure is a monotonic function of misregistration. This limited portion of the parameter space has been called the capture range. The correct registration transformation is, therefore, a local optimum value of the similarity measure within the capture range. As Studholme points out , this capture range is a function of the image content and field of view as well as of the similarity measure, so it cannot be determined a priori. An indication of the comparative sizes of the capture ranges of different similarity measures can be obtained by registering a set of test data with a large number of different starting estimates [23,24]. The similarity measure that converges to the correct registration solution for the largest range of starting estimates has the largest capture range for that data. An important consequence of the existence of a capture range is that the starting estimate of the registration transformation needs to be close to the correct transformation. For registration using mutual information, the capture range is often as large as a 30-mm translation or 30° rotation , which means that for the great majority of clinical images, registration can be fully automatic. For a minority of images, including situations in which the images being registered are acquired in different planes (e.g., one is axial and the other coronal), or where a small field of view in one modality reduces the capture range, a user might need to interactively align the data to provide the algorithm with an improved starting estimate. A further consequence of the capture range is that stochastic optimization algorithms such as simulated annealing and genetic algorithms are inappropriate for intermodality image registration using many voxel similarity measures. It is difficult to constrain these optimization algorithms to search just inside the capture range because neither the size of the capture range nor the position of the starting estimate within this range is known a priori.

5.4 Search Strategy

Calculus-based optimization techniques have already been described in the chapter entitled "within-Modality Registration using Intensity-Based Cost Functions". Most implementations of image registration by maximization of mutual information do not make use of derivatives. Instead they use standard numerical optimization techniques such as Powell and simplex, or simple multiresolution search. These algorithms need to find the optimum value of the similarity measure within the capture range, rather than the global optimum.

Powell's Direction Set Method

Powell's direction set method finds the optimum in an N-dimensional parameter space by a succession of one-dimensional optimizations. The advantage of this approach is that one-dimensional optimization is an easy problem to solve numerically. Each iteration involves carrying out N one-dimensional optimizations. For 3D rigid body registration, there are six degrees of freedom (three translations and three rotations) giving a six-dimensional parameter space. Determining transformations with more degrees of freedom involves a higher dimensional optimization space. For example, determining an affine transformation involves finding the optimum in a 12-dimensional parameter space. It is common to start by carrying out the one-dimensional optimizations along directions defined by the axes of the parameter space. For example, the method might start by searching in the lateral translation direction, then the anterio-posterio translation direction, then the cranial-caudal translation direction, then rotations about each of these axes in turn. This example gives only one of the 6! = 720 different orders in which the axes can be searched at the first iteration. The choice of ordering can influence the success of the optimization, and many researchers have adopted a heuristic approach to selecting the order of the search. Maes suggests  that for typical medical images that have higher in-plane resolution than through-plane resolution, it is best to search the in-plane degrees of freedom before the out-of-plane degrees of freedom. After an iteration, there is an estimate of the transformation solution Te that, in general, has a nonzero component for each degree of freedom. Subsequent iterations start their search at Te and can use either the same axes for the one-dimensional searches, or Te to define a new set of search axes. The algorithm finishes when an iteration fails to improve the similarity measure evaluation by more than some defined tolerance. A commonly used implementation is that provided in .

The speed of registration can be increased by applying the Powell algorithm at multiple resolutions. The resolution of the data is reduced by blurring, and then the blurred images are subsampled at intervals corresponding to the resolution. The algorithm is run at this resolution; then when it terminates, it is started again at a higher resolution (less blurring of the data), using the solution obtained at the lower resolution as the starting estimate. This increases the speed of the algorithm substantially both because the blurred and subsampled images have fewer voxels than the original, making the similarity measure quicker to calculate at each iteration, and also because fewer iterations are needed overall.

Simplex

The simplex optimization algorithm does not use a succession of one-dimensional optimizations. A simplex in N-dimen-sional space is a shape with N + 1 vertices, and is thus the simplest shape that requires that many dimensions to exist. For the purposes of N-dimensional optimization, the simplex can be thought of as a creature with N + 1 feet that tries to crawl into the correct minimum. One foot is the starting estimate, and the others are randomly generated within the search space, which should lie within the capture range. The similarity measure being optimized is evaluated for each foot, and the creature decides to move its worst foot to a better position. Because the simplex is the simplest creature that requires N dimensions to exist, all of the remaining feet lie in an N — 1 dimensional plane. The creature first tries to reflect its worst foot across this plane. If the new position it lands on is better than the one it came from, it tries going further in that direction. Alternatively, it tries to place its foot in an intermediate position. If none of these moves gives it a better solution, it tries moving all of its feet closer to the best one. It then once again selects its worst foot, and has another attempt. The simplex stops when all its feet are within a predefined range of each other. A simplex is likely to find the nearest local minimum in the search space. Like the Powell algorithm, the simplex algorithm can be applied at multiple resolutions.

Multiresolution Search

Because of the capture range issue describe above, Studholme has chosen to use a simple step-based optimization technique [21,23,24]. The starting estimate is assumed to lie within the capture range, then the similarity measure is evaluated at that starting estimate, and with a single increment and decrement in each parameter of the spatial transformation model. This gives 13 evaluations for a 3D rigid-body transformation. The best evaluation is chosen as the new estimate, and the algorithm iterates until no step can be found that improves the value of the measure. The translation step size St is chosen as approximately the resolution of the data, and the rotation step size 56 is chosen such that it corresponds to a translation of St at 60 mm from the center of rotation. This scheme is run starting at low resolution, and when the algorithm terminates at that resolution, the resolution is increased. At the highest resolution, the step size is further reduced to provide subvoxel registration solutions. This approach can be computationally expensive, but it is robust and easy to implement.