Fast Image Alignment Using Anytime Algorithms

1 Introduction

The need to align, or register, two images is one of the basic problems of computer vision. It can be defined as the task of finding the spatial mapping that places elements in one image into meaningful correspondence with elements in a second image. It is essential for data fusion tasks in medical imaging [Hajnal et al., 2001] and remote sensing (e.g. [ColeRhodes et al., 2003]). It is also widely applied in tracking and automatically mosaicking photographs [Szeliski, 2004].

One of the most straightforward and widely used approaches is referred to as direct image alignment. It works by defining a similarity measure, D, as a function of a reference image, and a template image warped by a transformation with some parameters, φ. The computation of D typically requires examining all pixels in each image. The alignment problem becomes that of finding the values of φ that Doina Precup School of Computer Science McGill University Montreal, Canada

∗Mr. Brooks was partially supported by an National Science and Engineering Research Council Post-Graduate Studies award.

dprecup@cs.mcgill.ca

maximize the chosen similarity measure. A number of optimization techniques for smooth functions, such as gradient descent, have been used for this problem, and provide good solutions on a wide range of image types. However, these approaches can be slow, which reduces their usefulness in time-sensitive applications such as real-time video registration (e.g. [Wildes et al., 2001]) and medical image registration during surgery(e.g. [Pennec et al., 2003]). It is possible to increase the speed of processing by using only a subset of the pixels to compute D, but this can easily lead to a reduction in accuracy and reliability. Determining the size of the subset to use is typically done in an ad-hoc fashion, or using heuristics which are applicable only to certain domains. Furthermore, since a different number of pixels may be needed at different stages in the optimization, a fixed subset is necessarily a compromise.

In this paper, we propose a deliberation control framework using anytime algorithms [Dean and Boddy, 1988; Horvitz, 1987] to arrive at a principled solution to the speed vs. accuracy trade-off in this problem. The first step is to learn the properties of the similarity measure under consideration, in terms of accuracy vs. computation time, by training off-line on image pairs for which the transformation parameters are known. Given a new pair of images to align, we then use this knowledge to determine the number of pixels that need to be considered at each step of the optimization. In this paper, we explore the effectiveness of this approach using two common similarity measures, mean squared difference and mutual information, and a gradient descent optimizer. We tested the algorithm on several types of images: images of everyday scenes, multi-modal medical images and earth observation data (i.e. Landsat and Radarsat images). In all cases, using a deliberation control approach is faster than computing the transformation using all the image data and gives more reliable results than simply performing the optimization using an arbitrary, fixed, percentage of the pixels.

The remainder of this paper is organized as follows. In Section 2 we review the image alignment problem and in Section 3 we review methods of deliberation control using anytime algorithms. The details of how deliberation control has been implemented in the context of image alignment are given in Section 4. Finally, Sections 5 and 6 describe our experimental setup, results and conclusions.

2 Image Alignment

Direct approaches to image alignment work by defining a function D that measures the similarity between a fixed reference image, R(x) and a new image T(W(x,φ)). Here, we consider the images R(x) and T(x) to be continuous functions of position, x, defined on some space of coordinates X. The coordinate space is warped by W(x,φ) which is a mapping from X to X, parameterized by φ (e.g., a translation or rotation). Since our images are actually sets of pixels located at integer coordinate positions, we use linear interpolation as needed to determine the values of T(W(x;φ)) when W(x;φ) is not an integer. The similarity measure, D, is thus a function of the transformation parameters φ, and the problem of image alignment becomes an optimization problem, which can be solved using many standard techniques, e.g., gradient descent, second-order methods, stochastic programming etc. All these techniques require the repeated calculation of D, and/or its gradient, ∇φD, at different points in the space of possible transformations. This is by far the most computationally intensive part of the process.

Recently, feature based approaches to alignment have seen considerable success, and can operate faster than direct approaches for many applications. Nevertheless, the direct approaches can yield higher overall accuracy and continue to be used as a final adjustment step [Szeliski, 2004]. Furthermore, it is difficult to match features reliably when the images in question are not of the same modality. Hence, direct approaches are the method of choice in applications such as medical imaging [Hajnal et al., 2001] and geomatics [ColeRhodes et al., 2003], where high precision is required and multimodal imagery is common.

It has long been known that numerical optimization approaches can be significantly accelerated by using only a subset of the pixels in the images to estimate the similarity function. For example, many implementations of mutual information (e.g. [Ibanez et al., 2005]) use a random subset of the image data. It has been suggested that it may be more efficient to include only pixels with high derivatives in the calculation [Szeliski, 2004]. However, the size of the subset to be used is fixed in an ad-hoc fashion. Unfortunately, any fixed size is usually too much for some regions of the parameter space and too little for others. Furthermore, different types of images behave differently, so it is difficult to come up with a subset size that is appropriate for all cases. Instead, we will use a deliberation control mechanism to choose how much computation to perform at each step of the optimization.

3 Deliberation Control with Anytime Algorithms

In many artificial intelligence tasks, e.g. planning, the quality of the solution obtained depends on the amount of time spent in computations. Hence, trade-offs are necessary between the cost of sub-optimal solutions and the cost of spending time doing further computation. This process, called deliberation control, has been investigated in the context of real-time artificial intelligence and a number of approaches have been proposed [Horvitz and Zilberstein, 2001]. Deliberation control methods rely on two key components: algorithms that support partial evaluation, and knowledge about how those algorithms perform after different amounts of computation.

A class of algorithms supporting partial evaluation are anytime algorithms [Horvitz, 1987; Dean and Boddy, 1988], which provide a solution when run for any length of time. The solution quality is guaranteed to improve with the amount of computation performed. Deliberation control strategies using anytime algorithms have been applied to both theoretical and practical problems including robot control [Vlassis et al., 2004], constraint satisfaction [Wah and Chen, 2000] and shape extraction in image processing [Kywe et al. , 2006].

To formulate an effective deliberation control strategy using anytime algorithms it is necessary to have meta-level knowledge of their performance as a function of the amount of computation performed [Dean and Boddy, 1988; Horvitz, 1987]. This knowledge is stored in a performance profile. Performance profiles may be based on theoretical knowledge of the algorithm, on empirical testing of its performance at different computation levels, or a combination of the two. In any case, the decision to continue computation will be based on an estimate of the accuracy of the current result, and an estimate of the potential improvement if the algorithm continues to run [Dean and Boddy, 1988; Horvitz, 1987; Larson and Sandholm, 2004].

The simplest type of performance profile is a static one, which predicts accuracy as a function of the amount of computation completed. However, for our problem of interest, this is equivalent to simply using a fixed, arbitrary percentage of the pixels. If feedback about the current run of the algorithm is available, it can be incorporated in a more sophisticated approach. A dynamic performance profile [Larson and Sandholm, 2004] uses feedback to estimate the accuracy as the algorithm progresses. It is described by two functions: ˆa = Pfwd( f , p), maps the percentage of computation completed, p, and the feedback parameter, f , to an expected accuracy, ˆ a. The other, ˆp = Prev( f ,a), maps a and f to the expected percentage of computation required, ˆp. Conceptually, these two functions are inverses. However, both have to be maintained in general, to facilitate the decision making. A controller can use these functions to gradually increase the amount of computation performed, until the estimated accuracy is adequate for the task.

4 Deliberation Control in Image Alignment

As mentioned in Section 2, the most computationally intensive part of image alignment is the repeated evaluation of the similarity measure D and its gradient ∇D. The optimization algorithm needs this information in order to take a step in parameter space towards the optimal setting. Note that the calculation only has to be accurate enough to ensure that the next step is correct; determining these values exactly is not necessary. Therefore, we propose to implement D and ∇D as anytime algorithms, and to learn performance profiles describing their accuracy at different levels of computation.

4.1 Anytime similarity measures

We implemented two popular image similarity measures as anytime algorithms. The resulting implementations have to be able to support partial evaluation, as well as to continue an interrupted calculation efficiently. To achieve this, we redefine each similarity measure, D(φ), as a function D(φ, p) of both the parameters, φ, and the percentage of pixels to be used, p. To avoid biasing the computation towards one area of the image, the pixels are processed in a random order.

Mean Squared Difference

The negative mean squared difference DMSD is one of the most common similarity measures [Szeliski, 2004], and is suitable for images of the same modality. In its original form, DMSD is simply an average of the negative squared differences in intensity between corresponding pixels (the negative here is simply to get a similarity, rather than a distance measure).

where N is the total number of pixels in the image. The gradient of this similarity measure is also easy to compute:

Note that both DMSD(φ, p) and ∇φDMSD(φ, p) can be updated incrementally in the usual fashion.

Mutual Information

The mutual information (MI) image similarity measure [Viola and Wells III, 1995] is useful for images of different modalities. Our anytime implementation is based on the efficient mutual information implementation proposed in [Th´evenaz and Unser, 2000], which relies on a B-spline windowed representation of the joint probability distribution of the intensity levels in the two images. Specifically, let bRk , where k = 1. . .K, be a set of K bins of width dR for the intensity values in the reference image starting at b R 0 =min x R(x). Similarly, let bTl be the bins for the intensity values in the template image, where l = 1. . .L, the bins begin at b T 0 =min x T(x) and have width dT . Then the unnormalized joint distribution is an array of size K×L. The entry Pkl is equal to the number of pixels in R for which the intensity falls in bin k and the intensity of the corresponding pixels in the transformed image falls in bin l:

where δ is equal to 1 if its two arguments are equal, and zero otherwise. Th´evenaz and Unser (2000) use instead a soft version to compute the entries in the table, based on B-spline Parzen windows. Similarly to their work, we will use:

where β0 and β3 are 0th and 3rd order B-spline Parzen windows respectively. The normalization factor of Pkl is:

This is true because B-splines satisfy the partition of unity constraint. The mutual information can then be computed using the usual formula:

Note that the second factor in the denominator above is just the intensity histogram of the original image, which is computed only once, before the optimization process.

Because a 3rd order B-spline is differentiable, the gradient of the joint histogram, ∇φPkl(φ, p), can be computed in the usual way as well, and stored in |φ| tables, each of dimension K×L. In [Th´evenaz and Unser, 2000] it is shown that when using the above formulation the derivative of DMI is:

where γ is the normalization factor γ = (dT⌊pN⌋)-1 . Thus, when the derivative is required, we can compute it from the joint probability distribution, and its gradient.

The algorithm maintains the unnormalized probability distribution, and its unnormalized partial derivatives, as described above, which can be easily updated when more pixels are added, because the table entries are simple sums. The distance measure and its gradient are then computed as needed from this table.

4.2 Performance profiles

Exploiting the partial evaluation possibilities of the anytime similarity measures requires dynamic performance profiles describing their expected accuracy at different computation levels and feedback values. The notion of accuracy must be developed in terms of the optimization method being used. In this paper, we used a simple steepest descent optimizer, which is described in [Ibanez et al., 2005] 1 . Carter (1993) has analyzed a similar class of optimizers and has proven their convergence using the following measure of relative error:

A dynamic performance profile requires a feedback parameter which indicates the progress of a particular calculation run and based on the above equation the gradient magnitude is an ideal candidate. Thus, we designed our performance profiles to be tables mapping computation level and gradient magnitude to accuracy. To construct them we sampled the gradient at different computation levels at many points in the transformation space. For each computation level, p, we grouped the gradient magnitudes into bins and computed the expected accuracy, ¯Ep, of the gradient for each bin as follows:

1 More specific information may be found in the class documentation http://www.itk.org/Doxygen/html/classitk 1 1RegularStepGradientDescentOptimizer.html

(where the φi are the sampled points in the transform space).

The optimizer uses this table to progressively increase the amount of computation performed until the estimated accuracy reaches its criterion for acceptability. The analysis in [Carter, 1993] indicates that significant computational gains can be made, with a small (ε ≈ 10%) reduction in accuracy. Therefore we choose to have the optimizer seek an expected accuracy of 90% for each gradient that it computes.

A simple example of how the table can be used to control computation is shown in Figure 1. This table can act as a performance profile where two functions, ˆa=Pfwd( f , p) and ˆp = Prev( f ,a) are implemented through simple lookup. For example, suppose an optimizer requires an accuracy of 98%. On an initial probe, f (which is in our case the current estimated performance) is 0.55. By examining the row applying to 0.6 and less, we get a prediction that with p = 8%, we will obtain the desired accuracy level (arrow 1). After performing 8% of the computation, however, suppose f is now 0.1. Thus the accuracy is only 93% (arrow 2) meaning more computation is required. The required p is now estimated as 16% (arrow 3), and so on. In our case, the parameter p will be the percentage of pixels processed in the image. The feedback parameter is the magnitude of the gradient.

5 Experiments

To test the anytime algorithm approach, a number of performance profiles were generated, and alignments were performed. Four classes of images (shown in Figure 2) with at least two image pairs each were used in the testing process. The first image class was typical digital photos (DP) (images a-d). Both a) and d) were self-aligned and a) was aligned affinely against several images of the same scene taken from different camera positions (images b,c) using both similarity measures. The second class of images (M1) (images e,f) were slices from T1-weighted magnetic resonance imaging (MRI) volumes which were self-registered using the DMSD measure. The third class of images (EO) are patches from georeferenced, orthorectified Landsat 7 and Radarsat imagery that were registered to each other using the mutual information measure (rows 3 and 4). The final class of images (M2) are slices from previously registered volumes in different medical imaging modalities, including T1 and T2 weighted MRI, and computed tomography (CT) (last row). These images were aligned to each other using the mutual information measure.

1 Images 2a, 2b, 2c, and 2d from K. Mikolajczyk http://www. inrialpes.fr/lear/people/Mikolajczyk/. Landsat and Radarsat images (2g,2h,2i,2j, 2k,2l) from Natural Resources Canada http://geogratis.gc.ca. Medical images (2e, 2f, 2m, 2n, and 2o) courtesy Montreal Neurological Institute

5.1 Generating performance profiles

Before performing any alignments, performance profiles for each combination of image class and similarity measure were generated off-line. These profiles were constructed using training image pairs of the same modalities as the ones to be aligned. This training data was separate from the testing data used later. For each training pair, the value of ∇D was computed at 4000 points in the transform space, at 12 different computation levels, and the profile was constructed using the method described in Section 4.2. The resulting performance profiles are shown in Figure 5.1. Note that each profile has a roughly similar shape. At large values of the feedback parameter, ||∇φD||, little computation is needed to get a good result. For smaller values, however, progressively more computation is needed. This agrees with our intuition; since the image noise level remains constant, small values are progressively harder to measure. It is the curved shape of these graphs that allows us to realize important performance gains. A constant fraction of pixels would inevitably be too many for some feedback values and too few for others.

Despite their basic similarity, however, there are important differences between the profiles. For example, note that for both similarity measures, the medical images require a much greater percentage of computation for an accurate result for a given gradient magnitude. The alignment tests, discussed below, reveal that these images require more computation to align successfully. Also note that the DMI performance profiles for both geographic data and medical images indicate that no less than 30% of the pixels will ever be used for these classes of images using this method.

Figure 3: Performance profiles: The percentage of computation required to achieve Ep = 90%.

5.2 Alignment tests

We implemented our approach by extending an existing, welltested, image alignment implementation. Specifically, we used an implementation similar to the example MultiResImageRegistration1 in [Ibanez et al., 2005], pp. 257-63, with appropriate transforms, and with extensions added to the optimizer and measures to support deliberation control. The experimental procedure was as follows. Three sets of 20 random starting positions were created. Each set was at a different effective distance from the identity transform in order to test the algorithms over the capture range of the optimizer. For each combination of image, similarity measure and algorithm, the true transform was composed with these starting positions, and the result was used to initialize the alignment process. Alignments were performed using the standard approach (labeled GD100 in the graphs), the standard approach using only a specified percentage of pixels (labeled GDx, if x% of the pixels is used) and our anytime approach. The computational efficiency was measured in terms of running time. Each reported time was obtained on a 1.9GHz AMD Athlon machine with 3GB of RAM.

To determine if any observed reductions in runtime came at the expense of performance, we also measured the quality and reliability of the algorithms. The quality of a result was measured by the root-mean-square (RMS) position error compared to the true transformation position. To calculate this, 50 randomly placed points in the unit square were scaled to the image extent. These were transformed both using the known true transformation and the computed transformation. The reported RMS is the square root of the mean of the squared coordinate differences for these points.

The reliability of each method was measured by the number of failed alignments. A particular run was considered to have failed when the registration converged to a transformation that lead to RMS pixel position errors greater than 5 pixels, or when the registration failed to converge to an answer at all. In practice, we found that this criterion was rarely ambiguous. The alignment would either yield results that were much better than this RMS value, or much worse. Failed runs were not calculated in the average times or registration error.

Since each algorithm being tested was run on the same set of test inputs, we used paired comparisons to test for significance. For each criterion, the null hypothesis, H0, was that the algorithms had the same performance as the original. The RMS errors and run times were compared using pairwise ttests, and the failure rates were compared using the McNemar test. To adjust for multiple comparisons, we used the Tukey method for the RMS and timing data, and the Bonferroni method for the failure rates [Howell, 2002]. All tests were performed at the 95% confidence level.

Adjusting for multiple comparisons can tend to accept H0 when it should be rejected, artificially bolstering our argument. Therefore, we report all three possible cases. Where all tests rejected H0 we conclude that performance differs. When H0 was rejected pairwise, but accepted when adjusted for multiple comparisons we consider the result ambiguous, and finally when H0 was accepted by all tests, we conclude that the data do not indicate a performance difference.

The experimental results are summarized in Table 1 a) (for MSD) and b) (for MI). The tables show the runtime of each method for each image pair, as well as for all runs combined. We also show the failure rate, and the RMS error with standard deviation. For the MSD measure, little can be concluded. All the algorithms under test show some improvement in speed, without significantly affecting failure rate or RMS error. For this measure, there is little to distinguish our method from simply reducing the number of pixels. However, the results for the MI measure highlight the advantages of our method. Reducing the number of pixels by a percentage frequently incurs a statistically significant loss of quality or reliability. The anytime method, however, delivers significantly faster times without sacrificing either the RMS error, or the failure rate. This is particularly apparent in the overall results (bottom row,Table 1-b).

An advantage of our approach is its adaptability. It performed the alignment of digital photos (DP) with the MI measure using a little more than 10% of the original running time, without changing the failure rate. In other cases, particularly the multimodal medical image registration (M2), more pixels seem to be inherently required to successfully align the images. Our method adapts to that requirement and maintains a low failure rate by increasing the computation performed.

6 Conclusions and future work

We proposed to use deliberation control methods in order to improve the efficiency of computer vision applications. We implemented such methods for the image alignment problem and showed a significant improvement in speed without degrading the quality of the results. Even when the performance gains are limited, a major advantage of this approach is that the number of pixels used is determined using a training process. Our results show that arbitrarily selecting a percentage of the image data to use for alignment will lead to very different results on different classes of images. Our method gives a principled way to determine how much of the image data needs to be processed to achieve reasonable results.

In the future, we plan to further investigate the problem of multi-modal image alignment using such algorithms. In this case, the amount of data is prohibitive for exact methods, especially if the data is volumetric, and if the registration has to be performed in real-time during surgery. We will investigate the use of more sophisticated methods for obtaining the performance profile and doing the deliberation. We will also look at other ways of highlighting the anytime aspect of similarity measures commonly used in computer vision.