Optimized automatic registration 3D

The Optimized Automatic Image Registration method in MIPAV determines a transformation that minimizes a cost function, which represents the quality of alignment between two images. The standard registration problem of computing a transformation that best aligns a reference image to a target image is formulated as a mathematical optimization problem of determining the global minimum of a cost function. The method evaluates the cost function at several different image resolutions, starting with the lowest resolution. Each step of increasing resolution utilizes the previously determined optimal transformation as a starting point and further refines its values. This method usually works well with the images of the same modality.

Background

The main approach for determining an optimal transformation is to calculate a cost function, determine how it changes as individual transformation parameters are varied, and find the parameters that minimize the value of the cost function. For instance, Figure 1 shows a plot of a sample cost function for a selection of transformation parameters; in each plot, a single parameter is varied while all other are kept constant.

Another fundamental step in this registration method is to resample and interpolate the reference and target images at several different resolutions (starting with a coarse voxel size and gradually increasing the resolution), while searching for the min cost function.

An advantage of this multi-resolution approach is that the initial optimization considerably reduces computational cost, since the number of sample points is substantially less. In addition, the overall alignment is easier to find at a coarse resolution. Improved accuracy of the transformation paremeters is achieved using the highest resolution image that has not been subsampled.

The resampling process may influence the computed value of the cost function; therefore, several different interpolation methods are incorporated with this technique. These methods are discussed in more detail later in Resampling Interpolation[[1]] Taking into account the resampling and interpolation criteria, the Optimized Automatic Registration method begins searching for an optimum transformation using the Powell method at a coarse resolution. Once a transformation has been determined that minimizes the cost function, the transformation is further refined by progressively increasing the resolution.

The algorithm offers various resampling interpolation methods combined with the Correlation Ratio, Least Squares, Normalized Cross Correlation, and Normalized Mutual Information as a cost function.

Figure 1. The plots of the Correlation Ratio cost function versus some of the individual parameter values. In each plot (A, B, C, and D), a single parameter is varied while the other are kept constant

Image types

You can apply Optimized Automatic Registration method to color, grayscale and black-and-white 2D and 3D images.

Outline of the method

Pre-optimization steps

The minimum resolution of the reference and target images is determined. For more information, refer to Blurring [[3]]
Both images are resampled and interpolated to create high resolution isotropic voxels. The following interpolation methods are available: trilinear (a default), B-spline 3-rd and 4-th order, 4, 5, 6 Lagrangian, and Windowed sinc. Refer to Interpolation methods used in MIPAV[[4]]
The center of mass (COM) for both images is then computed and one translation step is performed to align the COM. Refer to Calculating the center of mass [[5]]

Optimization steps

As a part of the transformation optimization process the images are subsampled by 8, 4, and 2 times.

levelEight optimization. Both images are subsampled and interpolated, so that each image is 8 times smaller. The parameters of the transformation are then systematically varied, where the cost function is evaluated for each setting. The parameters corresponding to the smallest value(s) of the cost function are then used as the initial transformation for the next level in the optimization.
levelFour optimization. The images are subsampled and interpolated, so that each image is 4 times smaller. And the transformation corresponding to the smallest value(s) of the cost function are determined and used as the initial transformation for the next level in the optimization.
levelTwo optimization. Similar processing as the previous two levels except the images are, first, subsampled and interpolated, so that each image is 2 times smaller.
levelOne optimization - 1mm resolution images are used in this step and the transformation is generalized to include 12 degrees of freedom.
Post-optimization

Using the same interpolation methods as described in Step 2 and the optimal transformation determined above, the method transforms the reference image into the same coordinate system of the target image.

Following Sections 1.2-1.5 describe in detail the algorithms and methods used to perform Steps 1-8 of the optimization routine.

Pre-Optimization

Blurring

The registration algorithm first, determines the minimum resolution for each dimension of the reference and target images. The reference image is blurred to the resolution of the target image if one of the resolutions in a dimension of the target image is 50% or larger than the corresponding resolution in the reference image. Likewise, the target image is blurred to the resolution of the reference image if one of the resolutions in a dimension of the reference image is 50% or larger than the corresponding resolution in the target image. It is possible (but not likely) that both the reference and target images will be blurred.

Resampling

Once this isotropic 1mm resolution images have been obtained, the subsampled versions are also created. The sub-sampling algorithm then simply keeps every n-th point (that is, 2, 4 or 8) on the lattice in each direction. Therefore, the new volume contains 1/n3 as many points as the original.

Resampling can be performed for x, y, and z axes for 3D images, or for only x and y for 2D images.

Calculating the center of mass

To calculate the center of mass (COM), the method uses the right hand convention in 3D coordinate systems (X, Y, Z). Thus, in the image space, the left hand corner of the image is set to (0,0,0). The x axis goes left to right, y axis goes top to bottom and z axis goes into the screen.

To calculate the COM, we, first, define the characteristics function of an object in an image to be as described in equation 7 [[7]]

Equation 7

Then, an area of an image can be calculated as

Equation 8

Figure 2. Calculating the center of mass (COM)

And the center of mass, denoted by (xCOM, yCOM) for 2D images (for 3D images, add the z coordinate also) is given by the 1-st moments of the object:

Equation 9

Process

The process used in this method performs the following:

Calculates the COM for both images as described above in [the center of mass]
Aligns the COM of both images, which is an initial transfomation that involves only translation;
For each resampled image (which is 8, 4, 2, and 1 times), determines the transform that minimizes the cost function.

levelEight optimization

levelEight optimization is primarily concerned with determining the rotational transformation that will lead to the globally optimal transformation, which corresponds to the global minimum of the cost function. levelEight optimization uses the lowest resolution images and coarse rotation angles evaluated at relatively large step sizes.

Optimization

This Optimized Automatic Registration method is based on the empirical observation that finding the correct orientation, or rotation, is the most difficult task in image registration, since the three rotation parameters are highly coupled and most of the erroneous registrations that have been examined have happened primarily due to an incorrect orientation. Therefore, the search concentrates on the rotational part of the transformation.

Specifically, the reference image is oriented, interpolated, and the cost function evaluated for coarse rotations angles (-30, -30, -30), (-30, -30, -15), ..., (-30, -30, 30), where the values represent the amount of rotation in degrees about the x, y, and z axes respectively. Since there are three rotation angles and each angle can contain five different values, there are 125 possible angle configurations.

For each angle configuration, a 4-DOF local optimization is also performed to find the optimal translation and (global) scale. By default, the initial values are set to

0 for translation
1 for scale.

The best 20% of the cost values and/or the corresponding angle configurations (candidate local minima) are stored in a vector of minima that is used as a starting point for a further optimization, which uses a smaller step size over a narrower range of angles. Refer to http://OptimizedAutomaticRegistration3D.html.

For each parameter setting corresponding to the top 20% of the cost function minima, the algorithm performs rotations over the fine grid (which is by default 15 degrees with 6 degrees step) and evaluates the cost function.

The method now finds all angle configurations that have the corresponding cost function values lower than their neighbors from the vector of candidate local minima.

For each of these sets of parameters a 7-DOF optimization is then performed, storing the results of the transformation and costs before and after optimization in a vector of minima. #See also Degrees of freedom [[8]]

Since the relative costs of each candidate solution may change at higher resolutions, where structures become more sharply defined, the top 20% of the cost function minima points are considered for the next higher resolution (4) stage, rather than just a single best solution.

By default, the algorithm uses Trilinear interpolation to transform images to this new orientation.

<div+>Figure 3. LevelEight optimization.

levelFour optimization

The levelFour optimization step uses the interpolated images subsampled by 4 to determine the transformation that minimizes the cost function starting with the transformations determined in levelEight.

The optimization determines a 7-DOF transformation that corresponds to the minimum value of the cost function. This transformation is then perturbed and the cost function is computed for these new settings. The perturbations correspond to six-degrees for each rotation parameter and a global scaling factor of 0.8, 0.9, 1.1, and 1.2. A type of transformation (translation and/or global scale) is specified by a user. It includes the following transformations: Affine-12 (a default), Rigid-6, Global Rescale-7, and Specific Rescale-9. Here, numbers indicate degree of freedom (DOF).

A vector of parameters and top 20% of the min cost function values are considered for the next step, which involves images subsampled by 2.

See also Cost functions used in MIPAV algorithms [[9]] and Degrees of freedom [[10]]

levelTwo

In levelTwo, the method uses the images subsampled and interpolated by 2. The cost function is evaluated using the transformation parameters corresponding to the top 20% minima obtained from the levelFour. It finds the best minimum, and then optimizes it with the 7-DOF and then 9-DOF transformation. The method then returns the best minimum after optimization.

Note: if a user has limited the degrees of freedom to six, as for the rigid transformation, there will only be one optimization run, with six degrees of freedom.

levelOne

This step uses the unsubsampled interpolated images and computes the value of the cost function for each parameter setting determined in the previous levelTwo optimization. The parameters corresponding to the minimum value of the cost function with this resolution are then further optimized for transformations with 7-, 9-, and 12- DOF, unless a rigid (6-DOF) transformation is desired.

What is saved. The method stores

The set of parameters or vector at which the minimum was reached;
The value of the cost function at that minimum;
And the matrix, which was the true input into the cost function and represents the transformation that gives the minimum cost of differences between the images.

Additional registration parameters

To better align the images, one can also include into calculation additional registration parameters, such as weighting coefficients and advanced parameters used for calculating the min cost function. For more information, refer to Applying the AlgorithmOAR 3D [11]

Resampling Interpolation

The Optimized Automatic Registration method resampes images using an interpolation scheme, where anisotropic voxels in the Z direction are resampled into isotropic 1mm cubic voxels. New voxel values are computed using a weighted combination of existing voxel values within a defined neighborhood of the new voxel location. The interpolation methods are available in this implementation of the registration technique, include the folowing:

For more information about the interpolation methods used in Optimized Automatic Registration, refer to Interpolation methods used in MIPAV.

Cost Functions

A similarity or cost function measures the similarity between two images. During the Optimized Automatic Registration the adjusted image V is transformed using the transformation functions described above. And the similarity S(U;Vt) between the reference image U and transformed image Vt is then calculated. The Optimized Automatic Registration method searches for the transformation that gives the smallest value of the cost function, which we assume is the transformation that also gives the best alignment.

The cost functions which are implemented in this method include:

Cost functions used in MIPAV algorithms include: Correlation ratio,Least Squares, Normalized Cross Correlation, and Normalized Mutual Information.

Powell algorithm

The Optimized Automatic Registration method uses the Powell algorithm to find the global minimum of the chosen cost function. For more information about the Powell algorithm, refer to http://math.fullerton.edu/mathews/n2003/PowellMethodMod.html.

Degrees of freedom

The number of independent pieces of information that go into the estimate of a parameter is called the degrees of freedom (DOF).

For more information refer to Cost functions used in MIPAV algorithms, Degrees of freedom.

Applying the AlgorithmOAR 3D

To run the algorithm, complete the following steps:

Open an input and reference image.
Call Algorithms > Registration > Optimized Automatic Registration.
The Optimized Automatic Registration dialog box appears.

In the dialog box,

Select the image in the Register to box;
Complete the rest of the dialog box options. Refer to Optimized Automatic Registration.
Press OK.

The algorithm begins to run and the Registering Images window appears with the status. When the algorithm stops running the registered image appears in a new image frame. Refer to Figure 4.

\|	Recommendation: First, use the default algorithm settings. And then, after completing the first registration run, modify them depending on the result you receive.

**Figure 4. The input image (A), reference image (B) and result image (C). Registration was performed using the default algorithm settings.**



A (top) B (middle) C (bottom)

Applying the Constrained AlgorithmOAR 3D

To run the algorithm, complete the following steps:

Open an input and reference image.
Call Algorithms > Registration > Constrained Optimized Automatic Registration.
The Optimized Automatic Registration dialog box appears.

In the dialog box,

- Select the image in the Register to box;
- Complete the rest of the dialog box options. See Figure 5
- Press the Advanced Settings button, and then complete the Advanced OAR settings dialog box. In this dialog box, you can limit the number of iterations and also the translation's range. Refer to Figure 6 for the dialog box options.

Press OK.

The algorithm begins to run and the Registering Images window appears with the status. When the algorithm stops running the registered image appears in a new image frame.

Recommendation: For the Constrained Optimized Automatic Registration algorithm, first, use the default algorithm settings. And then, after completing the first registration run, modify them (including Advanced settings) depending on the result you receive.

Optimized Automatic Registration dialog box options


Input options
Register	Select the reference image from the list. The list alphabetically displays all opened images.		\|
Degrees of freedom	Specifies the spatial model used to restrict the type of linear transformation being applied during registration. Select one of the following: rigid-6, global rescale-7, specific rescale-9, and affine-12.
Interpolation	Specifies the type of interpolation the algorithm will use when resampling. The list includes: · Trilinear · B-spline 3-rd order · B-spline 4-th order · Cubic Lagrangian · Quintic Lagrangian · Heptic Lagrangian · Windowed sinc See also <a href="InterplationMethods.html#wp999048">"Interpolation methods used in MIPAV"</a>.
Cost function<a name="wp1915683"></a>	Specify a cost function which will be used to evaluate registration. The list includes: Least squares, Correlation ratio, Normalized cross correlation, Normalized mutual information. For more information, refer to Cost functions used in MIPAV algorithms
Use the max of the min resolutions of the two datasets when resampling	If this option is chosen in the dialog box, the method uses the maximum resolution of the two datasets. It throws away some image information, but works faster. If this option is not chosen the algorithms uses the minimum of the resolutions when resampling the images. This can be slower, but does not lose the information.
Initiate registration process by applying Least Squares	This option provides a way to register images using the corresponding landmark points or VOIs placed in both images (you must put at least 4 landmarks). Note that this option works only for images that require rigid transformation e.g., rotation and/or translation, but it does not work on images that require scaling.
Rotations</a>
Use this option if there is a need to make manual adjustments to registration.
Apply same rotations to all dimensions	If checked, applies rotation that you specified to all dimensions. Otherwise, you should select the axis - X, Y, or Z - to which rotation should apply.
Coarse angle increment and Fine angle increment
options can be used for two pass registration. The first pass would apply a coarse-level rotation to quickly arrive at an approximate registration, and then the second pass would apply a finer-level rotation to continue from the point where the first pass registration has finished. By default, the coarse angle increment is set to 15 degrees and fine angle increment is set to 6 degrees.
Weighted Images	This gives weights for certain areas within the reference or transformed image. These areas can be specified using VOIs or reference files. Here, higher weights mean a greater impact in that area on the registration.
No weight	- the weighting factor is not used for registration.
Register area delineated by VOIs only			uses VOI to register a selected area.
Weight registration	If this option is checked, the following parameters become available: Choose ref. weight opens the dialog where you can select the reference file which contains the information about the weighting factor. Choose input weight opens the dialog where you can select the input file which contains the information about the weighting factor.
Output options
Display transformed image	allows to view the transformed image in a separate window. By default, this option is active.
Interpolation:	- select the calculation method, which is used to produce the output image from the array of transformation matrices and input image. You can choose among Trilinear, B-spline 3-rd order, B-spline 4-th order, Cubic lagrangian, Quintic lagrangian, Heptic lagrangian, or Windowed sinc. The default choice is trilinear.
Advanced options
Pressing the Advanced button opens the Advanced OAR Settings dialog box where you can specify the additional parameters which will be used to calculate the min cost function. Parameters are as follows:
Multiple of tolerance to bracket the minimum	Recommended values are 10-60.		\|
Number of iterations	Recommended are 1-5 iterations.
Number of minima from Level 8 to test at Level 4	By default this is the best 20%, but you can enter your own number here.
Subsample image for speed		subsamples the image.
Skip multilevel search. Assume images are close to alignment.
OK	Applies the algorithm according to the specifications in this dialog box.
Cancel	Disregards any changes that you made in this dialog box and closes it.
Help	Displays online help for this dialog box.
Figure 5. The Optimized Automatic Registration dialog box options

Advanced OAR settings for Constrained Optimized Automatic Registration 3D


Pressing the Advanced button in the main Optimized Automatic Image Registration 3D dialog box (see Figure 5) opens the Advanced OAR Settings dialog box where you can specify the additional parameters and (or) constrains (if Constrained AlgorithmOAR 3D was called) that will be used to calculate the min cost function. Parameters are as follows:
Multiple of tolerance to bracket the minimum		Recommended values are 10-60.
Number of iterations	Recommended are 1-5 iterations.
Number of minima from Level 8 to test at Level 4		By default, this is set to 5, but you can enter your own number here.
Translation search range	These options allow the user to specify the limits in translations range (in the X, Y, and Z dimensions), and therefore, constrain the optimized registration. The options appear active only if the Constrained Optimized Registration algorithm is called.
Limit translation range	If checked, this option limits the range of translations by the values (in mm) specified by the user.
Apply same translation to all dimensions	If checked, this option applies same limits to all dimensions.
Limits on X translation (or all translations)	Specify the translation limit in X direction, or in all directions if the Apply same translation to all dimensions option is checked.
Limits on Y translation	Specify the translation limit in Y direction.
Limits on Z translation	Specify the translation limit in Z direction.
Subsample image for speed		subsamples the image.
Skip multilevel search. Assume images are close to alignment.
Initialize registration process by aligning COG's	If checked, the algorithm, first, aligns centers of gravity (COG) or centers of mass (COM) of both images, and then proceed with registration. See also Calculating the center of mass
OK	Applies the algorithm according to the specifications in this dialog box.
Cancel	Disregards any changes that you made in this dialog box and closes it.
Help	Displays online help for this dialog box.
Figure 6. The Advanced Optimized Automatic Registration (constrained) dialog box options

References

Bjórn Hamre "Three-dimensional image registration of magnetic resonance (MRI) head volumes" Section for Medical Image Analysis and Informatics Department of Physiology & Department of Informatics University of Bergen, Norway.

Bourke Paul "Interpolation methods" ~pbourke/other/interpolation/index.html http://local.wasp.uwa.edu.au/ ~pbourke/other/interpolation/index.html

B-splines: node4.html http://www.cl.cam.ac.uk/teaching/1999/AGraphHCI/SMAG/ node4.html

Derek L.G.Hill and Philipe Batchelor, Medical Image Registration, Chapter 3: "Registration Methodology: Concepts and Algorithms" edited by Joseph V. Hajnal, Derek L.G.Hill, and David J. Hawkes,CRC Press, 2001, pp. 40-70.

Roger P. Woods, Chapter 33 "Within-Modality Registration Using Intensity-Based Cost Functions" in Handbook of Medical Image Processing and Analysis, Editor, Isaac N. Bankman, Academic Press, 2000, pp. 529-553.

FLIRT homepage: http://www.fmrib.ox.ac.uk/fsl/flirt/

Jenkinson, M. and Smith, S. (2001a) "A global optimisation method for robust affine registration of brain images". Medical Image Analysis, 5(2):143-156, http://www.fmrib.ox.ac.uk/fsl/flirt

Jenkinson Mark and Stephen Smith "Optimisation in Robust Linear Registration of Brain Images" FMRIB Technical Report TR00MJ2, http://www.fmrib.ox.ac.uk/fsl/flirt

Josien P. W. Pluim, J. B. Antoine Maintz and Max A. Viergever, "Mutual information based registration of medical images: a survey", IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. XX, NO. Y, MONTH 2003.

Powell method to find the global minimum: mathews/n2003/PowellMethodMod.html http://math.fullerton.edu/ mathews/n2003/PowellMethodMod.html

Stephen M. Smith , BET: "Brain Extraction Tool", FMRIB technical Report TR00SMS2b, FMRIB (Oxford Centre for Functional Magnetic Resonance Imaging of the Brain), Department of Clinical Neurology, Oxford University, John Radcliffe Hospital, Headington, Oxford OX3 9DU, UK.

Thomas M. Lehmann, Claudia Gonner, and Klaus Spitzer. Survey: Interpolation Methods in Medical Image Processing. IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 11, NOVEMBER 1999 1049.