2D and 3D Multimodal Image Registration

by

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.

- It assumes that the images being registered are similar in intensity up to a geometric transformation (i.e. scaling, rotation and translation).
- It assumes that image intensities between the two images will be maximally correlated for the "correct" transformation.

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:

- Minimize the squared differences between corresponding values in tape 1 and tape 2.
- For each value in tape 1, minimize the variance of corresponding values of tape 2.

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 **Y _{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.

**Y _{T}** denotes

**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 **Y _{T}**.

We must further define the mean pixel value for the image **Y _{T}**as:

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 **Y _{T }**. To do
this, we must create a mapping which allows us to associate a given pixel
in image

This joint probability is defined in "intensity space" where
there are two orthogonal axes which represent the intensities in the images
**X _{L}**and

This is essentially a "2D histogram" of the image pair.
We can now marginalize this joint distribution by summing over all possible
values of **X _{L}**. To calculate the marginal variance, we
consider each intensity of

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 **X _{L} **is calculated using

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 **X _{L} **is calculated using

**Experiment 1**

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.

**Experiment 2**

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.

**Unperturbed 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.

- The minimization algorithm is converging to a local minima.
- The images obtained are not actually in correct registration.
- The correlation ratio may assume a maximal value for a transformation
which is not
*exactly*equal to the position of correct registration.

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.

**Acknowledgments**

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**

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.