2D and 3D Multimodal Image Registration
Maximization of the Correlation Ratio
EN292 - Medical Imaging Final Project
by: Andrew Willis
May 17, 1999
Recent advances in computing power have made it possible to use data provided from medical images in new and revolutionary ways. A crucial aspect of these developments hinge upon coalescing image data provided from different modalities. This is done by image registration. Given two images representing the same or similar structures, we want to automatically determine the transformation which will allow these structures to be superimposed. The complications of this problem has sparked a line research on the topic of "image registration."
I have chosen to implement a registration method recently suggested by a group of researchers from INRIA, a French research agency. Alexis Roche, Gregoire Malandain, Xavier Pennec, and Nicholas Ayache have proposed a method of registration based on the maximization of a function called the "Correlation Ratio." I have implemented this method for the registration of 2D images and analyzed it's performance in applications between several modalities.
This method of image registration makes some assumptions regarding the two images to be registered.
The foundation of the correlation ratio is based upon viewing image intensities as random variables. By using statistical measures of the image intensities are results are not based directly upon the image intensity values. Rather, we are using the mean and variance of the pixel intensities as a measure of comparison.
Consider the following problem:
We are given the ticker tapes of numbers and we want move these ticker tapes laterally to get the best match for the number 2. We will consider two approaches:
For value 2 the squared difference is 34, the variance for value 2 is approximately 6.
For value 2 the squared difference 75, the variance for value 2 is 0.
Considering only value 2 there are two orientations where the variance for value 2 is 0. Consequently, we consider the variance for all values in the top tape (1, 2, and 3). We assume that minimizing this total variance will result in the best possible alignment of the two tapes.
This is the underlying assumption of the correlation ratio and forms the entire framework for correlation ratio based registration.
The correlation ratio is a "similarity measure" between images X and YT for a given transformation T.
To thoroughly understand the mechanisms of this registration functional we will consider each of the constituent variables and equations.
X denotes the template image. We fix the position of this image.
YT denotes floating image. This describes the image resulting from putting image Y through the transformation T.
T denotes a rigid transformation which involves both rotation and translation.
In layman's terms:
This equation is an expression of the statistical variation of intensities in the image Y after a transformation T (see specifics).
N denotes the total number of intensity samples in image YT.
We must further define the mean pixel value for the image YTas:
In layman's terms:
These equations completely describe how to compute the denominator of the correlation coefficient.
Let us now concentrate on the numerator:
Calculation of this value requires the definition of the joint probability of intensities in images X and YT . To do this, we must create a mapping which allows us to associate a given pixel in image X with a pixel in YT. This is done by subjecting image Y to a given transformation T and then linearly interpolating the image X to align it with the image YT. Linear interpolation of image X creates a slightly different image we will denote XL. Once the two images are aligned we can pair the intensity of corresponding pixels in the two images together. These pairs are then used to define the joint probability distribution of the intensities in images XL and YT (given a rigid transformation):
This joint probability is defined in "intensity space" where there are two orthogonal axes which represent the intensities in the images XLand YT.
This is essentially a "2D histogram" of the image pair. We can now marginalize this joint distribution by summing over all possible values of XL. To calculate the marginal variance, we consider each intensity of YT and calculate the mean and the variance of the paired intensities in image XL (see specifics). We then calculate the average value of the variance over all of the intensities which gives us the value of the numerator.
Now that we have totally defined how to calculate the correlation ratio given two images X, Y and a transformation T. We would like to find the value of T which maximizes the value of the correlation ratio:
A multitude of methods are available for maximizing a function.
In practice, we are actually trying to minimize the average variance functional which forms the numerator of the correlation ratio equation originally stated. The minimization method implemented for this project was the downhill simplex method in multidimensions.
This registration method may be applied in 2D and 3D with little theoretical modification. This section will discuss briefly the specifics involved with implementation in 2D and 3D.
In 2D registration we are given two 2D images which we would like to bring into alignment. For this project, I assume that the scaling factors to convert the images to the same pixel size is given. As a pre-processing step, the larger of the two images is subsampled to have the same pixel size as the smaller image. This creates two images of equal dimensions X and Y.
Now, image Y is subjected to transformation. In this case, the rigid transformation is given by:
In this case our transformation vector T is made up of 3 unknown parameters which we are attempting to find:
We then repeatedly calculate the correlation ratio for different transformations until we find the optimal transformation (i.e. the transformation that maximizes the correlation ratio).
The interpolated image XL is calculated using bilinear interpolation.
In 3D registration we are given two 3D images which we would like to bring into alignment. For this project, I assume that the scaling factors to convert the images to the same pixel size is given. As a pre-processing step, the larger of the two images is subsampled to have the same pixel size as the smaller image. This creates two images of equal dimensions X and Y.
Now, image Y is subjected to transformation. In this case, the rigid transformation is given by the following 4 equations:
Rotation about Z axis
Rotation about X axis
Rotation about Y axis
Translation from origin
In this case our transformation vector T is made up of 6 unknown parameters which we are attempting to find:
We then repeatedly calculate the correlation ratio for different transformations until we find the optimal transformation (i.e. the transformation which maximizes the correlation ratio).
The interpolated image XL is calculated using trilinear interpolation.
To check the registration algorithm I replicated a synthetic experiment performed by the INRIA authors. In this experiment a 2D image was subjected to the following intensity transformation. For each pixel in image X the corresponding pixel in image Y was calculated using the following formula:
Let us watch the registration procedure for this experiment. In this experiment, transformation parameters are randomly generated.
This experiment consistently resulted in very accurate subpixel registration.
A second experiment was performed where each of the pixel values in a T1 MRI image were inverted:
Y(x,y) = | X(x,y)-255 |
We then register this image to itself. Let us watch the registration procedure for this experiment. In this experiment, transformation parameters are randomly generated.
This experiment consistently resulted in very accurate subpixel registration.
An experiment was performed where MRI T1 and T2 images were obtained from the UCLA image database. Corresponding slices of T1 and T2 image data were taken from the database. The images were then registered without any transformation of the images.
Let us watch the registration procedure.
The result always has a negative translation in the Y direction. This is easily explained by closely looking at the images. One can see the obvious difference in intensities at the bottom of the images (back of the skull). This leads to a basic question. Are the images correctly registered or is the correlation ratio assuming a maximal value at a position other than the optimal registration position?
Look at the original images and see if you can discern a verticle displacement. I believe that one may exist.
T1 Image - Slice 35 T2 Image - Slice 35
The T1 and T2 MRI images obtained from UCLA had been registered to one another via a landmark based method. Details of this are outlined at the UCLA image database.
MRI T1 and T2 images were obtained from the UCLA image database. Corresponding slices of T1 and T2 image data were taken from the database. The T2 image was then subjected to a rigid transformation.
Let us watch the registration procedure for T1 - T2 Registration. In this experiment, transformation parameters are randomly generated.
This experiment consistently put the images a close proximity to register. However, I found that the registration did not exactly correspond to the original transformation. This may be explained by my previous experiment.
MRI T1 and PET images were obtained from the UCLA image database. Corresponding slices of T1 and PET image data were taken from the database. The PET image was then subjected to a rigid transformation.
Let us watch the registration procedure for T1 - PET Registration. In this experiment, transformation parameters are randomly generated.
This demonstrated that registration between T1 and PET images is difficult using 2D images. Seldom are the images properly aligned. In this case, I believe that additional data and initialization closer to the correct positioning via HCI (human computer interaction) would yield better results.
I found that for generalized multimodal image registration, the correlation ratio registration requires more data than that provided by 2D slices. As shown in the T1 - PET registration examples. However, in images which have some common structure evident, registration can be done consistently and reliably. This was demonstrated by the T1-T2 registration experiment.
It is also important to note that this algorithm, as with most all automatic registration algorithms, depends upon the minimization of an energy functional. Consequently, it is subject to the possibility of falling into local minima which give erroneous registrations. The sensitivity of the algorithm to the effects of local minima may be reduced by providing some rough alignment information which will bound the parameter space (i.e. roughly align the images using externally provided information from an HCI).
Analysis of the unperturbed experiment suggests to me that it is possible that the original monkey data which I used may not be exactly aligned as originally thought. Also, if we analyze results obtained from T1 - T2 registration, we can see an important issue. The solution is seldom exactly equal to the initial transformation (within subpixel accuracy). There are three possible sources of error.
My experiments indicate that the minimization algorithm is performing acceptably. In the unperturbed experiment, I believe that the correlation ratio may have correctly registered the MRI images where the initial landmark based method used at UCLA may not have been completely accurate.
Naturally, the abundance of data available in a 3D data set will enhance this registration technique and should provide more accurate and robust results.
Comparison has been done with a landmark based registration and differences have been identified. A stereotactically aligned dataset may resolve the ambiguity as to which method is more accurate to the actual optimal registration position.
Extension of this algorithm to 3D is directly motivated by the conclusions of this project.
My minimization algorithm utilizes the downhill simplex method. There are other methods which may perform better in locating a global minima.
I would like to thank the people at the UCLA Laboratory for Neuro-Imaging whom have generously allowed public access to their image database.
Source code was developed for implementation over the web. I have written some supporting code which is particularly useful if one wishes to use Java for any serious image processing.