Reducing interpolation artifacts for mutual information based image registration.

Medical image registration methods which use mutual information as similarity measure have been improved in recent decades. Mutual Information is a basic concept of Information theory which indicates the dependency of two random variables (or two images). In order to evaluate the mutual information of two images their joint probability distribution is required. Several interpolation methods, such as Partial Volume (PV) and bilinear, are used to estimate joint probability distribution. Both of these two methods yield some artifacts on mutual information function. Partial Volume-Hanning window (PVH) and Generalized Partial Volume (GPV) methods are introduced to remove such artifacts. In this paper we show that the acceptable performance of these methods is not due to their kernel function. It's because of the number of pixels which incorporate in interpolation. Since using more pixels requires more complex and time consuming interpolation process, we propose a new interpolation method which uses only four pixels (the same as PV and bilinear interpolations) and removes most of the artifacts. Experimental results of the registration of Computed Tomography (CT) images show superiority of the proposed scheme.

1. INTRODUCTION

Image registration is the alignment of two or more different images of a view or an object such that their similar pixels coincide to each other. These images might be acquired by several sensors or by a single sensor in different viewpoints, moments or other different conditions. The goal of image registration is to find a transform function by which the floating image is aligned to the reference image. Every image registration algorithms require three basic ingredients: 1) a spatial transformation model which determines the set of possible solutions, 2) an objective similarity measure which estimates the quality of each potential solution, 3) an optimization algorithm which looks for the best solution[12].

There are numerous medical imaging modalities that show the anatomy or morphology. Especially, CT and magnetic resonance (MR) imaging of the head for diagnosis and surgical issues, provide significant information for surgeons. To obtain more complete information about the patient, monitoring tumor growth and comparison of the patient's data with anatomical information, registration between these modalities is necessary[23].

Multimodality registration is of great importance in many medical applications. Over the years, many different methods are advised for this task. A complete survey of these methods and classification of approaches can be found in[45].

Image registration methods are classified to 1) intensity and 2) feature based methods. Feature based methods use some features (edge, surface and line etc), which requires segmentation and feature selection. Intensity based methods, which are more popular in medical image registration[6], employs the intensity of pixels and no feature is needed.

Mutual information is an information-theoretic concept which has been widely used as the similarity measure in intensity based methods. It shows high accuracy and robustness with intensity variation and noise[78]. Mutual information is also employed in other image processing applications such as template matching[9] and object tracking[10].

The most important and critical step in calculating the mutual information of two random variables (or two images) is to find their joint probability distribution. Since the acquired images are actually samples of continuous two-dimensional signals, calculation of their joint probability distribution cannot be perfect. Thus an interpolation method should be employed for estimating the joint probability distribution which influences the quality of image registration. Therefore, selection of a reliable interpolation method for estimating the joint probability distribution is of great importance.

The joint probability distribution can be derived by using PV or bilinear interpolation methods[11]. When the underlying images are acquired by a single sensor and hence the size of their pixels is equal, mutual information (as a function of translation) shows some local optima at integer values[12]. These artifacts may cause the optimization algorithm to converge to wrong solutions. This will degrade the accuracy of the registration.

In order to remove such artifacts, HPV and GPV interpolation methods are proposed[513]. With these methods interpolation procedure is performed by using Hanning windowed sinc and Spline as kernel functions instead of boxcar function used in PV interpolation. In[5], it is shown that the performance of these methods is better than that of PV. But this comparison is questionable, since in HPV and GPV 16 pixels are involved in the interpolation while for PV only 4 pixels are involved. As a result, one cannot conclude that the superiority of HPV and GPV is due to their appropriate kernel functions. In contrast, in this paper we show that the importance of the number of involved pixels in the performance of registration is more than the selected kernel function. Also we propose a new interpolation method which uses only four pixels with better performance with respect to PV and bilinear interpolations. Since the resulted mutual information function is smoother, the registration is more accurate.

The paper is organized as follows. In Section 2 the problem of generated artifacts in mutual information function is stated. In Section 3 different interpolation methods are explained. The main idea of the paper is also discussed in Section 3. Finally, in Section 4 performance of PV, Bilinear and proposed method are compared.

2. ARTIFACTS IN MUTUAL INFORMATION FUNCTION

Mutual information is one of the basic concepts of Information theory which indicates the dependency of two random variables. Initially, it was used in medical image registration by Viola and Wells in 1995[14].

It is well-known that two random variables A and B with marginal probability distributions p(a) and p(b) and joint probability distribution p(a, b) are independent iff p(a, b) = p(a).p(b). Mutual information I(A, B) measures the dependency of A and B through the Kullback-Leibler distance of the distributions p(a, b) and p(a).p(b)[15], i.e.

Probability distributions of images are easily derived by normalizing the correspondent histograms. Mutual information takes its maximum value when the underlying random variables (or images) are absolutely dependent (one of them is a function of the other one). In such case, the joint probability matrix (or joint histogram matrix) is diagonal. Reversely, if the images are independent, then the mutual information takes its minimum value, i.e. zero.

In order to align two images, a transformation function T should be applied to the floating one. Transforming a grid point of floating image usually gives a point in the reference image which does not coincide with a grid point. Thus an interpolation method is required to estimate the joint histogram of the images.

The prevalent interpolation method in image processing is bilinear interpolation which uses four nearest pixels of the transformed coordinate in the reference image. Joint histogram is easily obtained by a simple counting procedure[11]. With this method the resulted mutual information function is not smooth and the optimization algorithm may converge to a wrong solution. The situation is more critical when the images are noisy.

In order to overcome this problem, PV interpolation method was proposed. But it has been shown in[12] that if the underlying images have the same pixel size, PV and bilinear interpolation yield some local optima (artifacts) on the surface of mutual information function. A typical pattern of mutual information function in terms of displacement (translation) in horizontal and vertical directions is shown in Fig. 1. It can be seen that the surface has some local minimums (res. maximums) for bilinear (res. PV) interpolation method.

Figure 1

Typical patterns of local extreme resulting from two interpolation methods. (a) Mutual information value for x and y translation by PV interpolation; (b) Mutual information value for x and y translation by bilinear interpolation

3. INTERPOLATION METHODS AND PROPOSED ALGORITHM

In this section, first we describe some interpolation methods such as bilinear, PV and its modifications. We investigate the reason behind the acceptable performance of GPV and HPV methods. Also, an interpolation scheme is proposed.

A. Bilinear Interpolation

As mentioned in the previous section, transforming coordinate of a pixel of floating image F gives a new coordinate in the reference image R, which may not be a grid point. Hence the intensity of such point must be somehow interpolated. In bilinear interpolation, the intensity is estimated by computing a weighted average of the intensity of four nearest pixels in the reference image.

If the point p with coordinate s in the floating image is transformed (by transformation T) to point q with coordinate T in the reference image (see Fig. 2 for more details), the intensity of the reference image in coordinate T is given by

Figure 2

shame of transformed pixel from floating image which conceded to a point with coordinate TS and four nearest pixels in reference image.

where n, i=1,2,3,4 denote the coordinate of four nearest grid points to point q in the reference image. Note that we have

For each point p of the floating image the command h(R(T,F(s))+ = 1 is used to estimate the joint histogram of floating and reference images. Since bilinear interpolation does not show acceptable robustness to noise[6], other methods such as PV are more popular.

B. PV, GPV, HPV and the proposed interpolation methods

PV method is not actually an interpolation scheme for estimating the intensity in a non-grid point. It is only designed to estimate the joint histogram of two images. Instead of interpolating the intensity values in R, the contribution of the intensity of the point p to the joint histogram is distributed over the intensity of all four nearest neighbors of q. Updating the joint histogram of F and R is performed by h(R(n),F(s))+ = w where

w = ƒ(d1)ƒ(d2) , i = 1,2,3,4      (3)

and n, i=1,,2,3,4 denotes the four nearest grid points to point q in the reference image. Also d1, d2 stands for the horizontal and vertical distance of q from n.

With the PV method, the kernel function is defined by:

ƒ(x) = 1 – x , 0

Which is known as Boxcar function. In[9] Chen and Varshney proposed GPV interpolation method. They replaced Boxcar function with a third order Spline, i.e.,

In this method sixteen pixels are used to calculate w 's and update the joint histogram. As a result, for each pixel of floating image, sixteen bins are updated in the joint histogram. Note that the support of this function is the interval[0,2].

In[8], HPV interpolation method was proposed by replacing the boxcar function with Hanning windowed sinc function

Similar to GPV method, sixteen pixels are used to estimate the joint histogram. The kernel functions f, f, f are shown in Fig. 3.

Figure 3

kernel functions. (a)-Boxcar. (b)- Third order Spline. (c)- Hanning windowed Sinc. (d)-, (e)- and (f)-

It can be seen from that their behavior are almost the same. Thus we can expect to achieve the same performance with other kernels with such behavior. For instance, we propose the functions f1 and f2 (shown in Fig. 3)

Using these kernel functions for updating joint histogram, the mutual information function is plotted in Fig. 4 as a function of translation and rotation. All methods which use sixteen pixels (i.e. GPV, HPV, f1 and f2) for interpolation have the same pattern. In particular, no local optimum is seen on the surface of the objective functions. But bilinear and PV interpolations that use only four pixels, are deeply affected by the artifacts.

Figure 4

Typical patterns of mutual information versus rotation and translation using bilinear, PV, HPV, GPV, f1 and f2 (from first to sixth rows respectively).

Here we use another kernel function

f3(x) = 1 – sin(πx2 / 2)

with four pixels. Although its behavior is similar to four mentioned ones (see Fig. 3.d), its associated objective function is not smooth (like PV interpolation). Therefore we can conclude that this is not the special kernel function which removes the artifacts in GPV and HPV methods. In fact, the number of involved pixels (16) is much more significant.

On the other hand, using the proposed f2 function (shown in Fig. 5) shows smoother surface in noisy images. Thus it can be claimed that if sixteen pixels are used for interpolation, there are kernel functions which removes the artifacts even better than HPV and GPV methods.

Figure 5

Patterns of mutual information versus translation using three different methods for noisy image.(a)-HPV. (b)-GPV. (c)-f2

Since using sixteen pixels for interpolating is a time consuming task, we are motivated to propose a new scheme which performs better than PV and bilinear interpolation methods. The distance of transformed pixel in the floating image from its four nearest pixels in the reference image is directly used in PV and bilinear interpolation methods which yields to some local optima on the surface of the mutual information function.

We try to eliminate such local optima by using the intensity of four nearest pixels besides their distances from the transformed pixel and changing the kernel function. We use f(x) = exp((1 – x2)/2) as the kernel function and multiply it by the intensity of the four nearest neighbors.

We define

Where R(n) is the intensity at point n and d1, d2 denote horizontal and vertical distance of the transformed point from n1. Similar to PV relation, (3) is used to update the histogram. Note that if we use only the intensity of the pixels to calculate the joint histogram (the kernel function in equation (9) is ignored), local optima will disappear but this type of interpolation, has no strong peak in the ground truth (see Fig. 7.a).

Figure 7

Pattern of mutual information in terms of translation. (a)- registration between images Fig.6.a and Fig.6.b by using only intensity of pixel for interpolation. (b)- registration between images Fig.6.a and Fig.6.b by proposed interpolation method. (c)- registration between images Fig.6.a and Fig.6.c by proposed interpolation method.

Figure 6

(a)- Original image (b)- original image translated by 4 pixels (c) original image translated by 4 pixels which its intensity is decreased.

It may seem that direct use of intensity in equation (9), makes the similarity measure sensitive to intensity variation. Certainly the joint histogram of the images (and hence their joint probability distribution and mutual information) will be changed when the intensity of the pixels changes. This happens for all known interpolation methods. The key point is that the mutual information is robust enough (to intensity variations) so that the accuracy of registration is not affected (see Fig. 7.b and Fig. 7.c).

Fig. 8 shows mutual information function in terms of rotation and horizontal translation by PV, bilinear and proposed interpolation methods. It can be seen that the local optima are removed with the proposed method. However, the surface of objective function is not smooth in PV and bilinear interpolation. As a result, the accuracy of registration is better with the proposed method.

Figure 8

Patterns of mutual information versus translation and rotation using three different methods. (a)-bilinear. (b)-PV. (c)-proposed method

4. EXPERIMENTAL RESULTS

For evaluating and comparing the performance of our method with PV and bilinear methods, rigid body transformation model was selected. With this assumption only rotation and translation are permitted. Thus for 2D images, three parameters (rotation around z axis and translation in x and y axis) are considered.

Note that comparison between as all of them use 4 pixels in interpolation. We used fifth slice of CT images (512×512) of five patients from dataset of the project “Retrospective Image RegistrationEvaluatio”. Applying the transform function T(θ,t,t), each images was rotated θ degree around z axis and translated in x and y directions by t and t pixels.

Thus in this experiment mutual information metric is a function in terms of three variables (θ,t,t). Zero-mean Gaussian noise with variance 0.1 was also added. The number of histogram bins in all images was 128. Mutual information was used as the similarity measure and Powell-Brent method employed for optimization. Powell-Brent method is one of the popular methods[16]. With this method, in each step of the optimization process mutual information function is optimized in one direction of θ, t or t and the resulted point in each step is used as the starting point for the next step. In the next iteration of optimization the resulted points in the previous iteration are used. In our simulation, point (0, 0, 0) was considered as the starting point and variation of parameters (θ,t,t) was limited to the interval [-10,10] . The result of this simulation is shown in Table 1.

TABLE 1

THE RESULT OF REGISTRATION USING THREE DIFFERENT INTERPOLATION METHODS

It can be seen that both of the proposed methods and PV interpolation give estimation of transformation function T better than that of the bilinear interpolation. But when the images are not rotated (or θ is zero), the proposed method has high accuracy and fewer errors. This is because the PV interpolation artifacts are the direct result of using the simple boxcar kernel function and when θ=0 local optima will appear for integer pixel shifts. Such artifacts are not so evident when θ is not zero[17]. As seen in Fig. 6, local optima are in integer translations and the resulted estimations for patients 3, 4 and 5 (θ is zero for these patients) which are given in Table 1, confirms this.

To clarify the superiority of the proposed method we set θ=0 in another experiment. Original images were added by zero-mean Gaussian noise with variance 0.01. Optimization was performed to drive optimal value of t and t. The results of this simulation are given in Table 2.

TABLE 2

THE RESULT OF REGISTRATION USING PV AND PROPOSED INTERPOLATION METHODS

It can be seen that with PV interpolation method, the optimization algorithm converges to one of the local optima (in integer translation) which are shown in Fig. 6. However, due to removing artifacts, the proposed method shows higher accuracy.

5. CONCLUSION

In this paper the intensity-based image registration methods which use mutual information as similarity measure and the effect of artifacts resulted by PV and bilinear interpolation methods was studied. Also we considered GPV and HPV interpolation methods that use sixteen pixels for interpolating. By introducing two new functions f1 and f2 as kernel functions, we showed that the acceptable performance of GPV and HPV is because of using sixteen pixels not using special kernel functions. Since methods which use sixteen pixels for interpolation are time consuming, we proposed a new interpolation method that use four pixels to interpolate images and removes local optima (artifacts) from the objective function.

Hossein soleimani received the B.sc. and M.Sc. degree from Shahid Rajayi Teacher training University, Tehran, Iran, in 2008 and Isfahan University of Technology in 2011 respectively, both in electrical Engineering. During his four years at Shahid Rajayi Teacher Training University his research was in neural network systems and image processing. His interests in Isfahan University of Technology include image and signal processing and medical imaging.

Mohammadali Khosravifard received the B.Sc., M.Sc., and Ph.D.degrees in electrical engineering from Shiraz University, Sharif University of Technology, and Isfahan University of Technology, Iran, in 1996, 1998, and2004, respectively. Since 2004, he has been an Assistant Professor with the Department of Electrical and Computer Engineering, Isfahan University of Technology. His main research interests include information theory and related topics in image processing.