EJBI - European Journal for Biomedical Informatics :: Low-dimensional Multimodal Deformable Registration of MRI Brain Images in Stereotaxic Space

One of widely applied methods in computational neuroanatomy is a voxel-based morphometry (VBM), which has recently become a subject of discussion [1], [2]. It interrogates anatomical MRI scans on voxel by voxel basis, in order to demarcate regions with significant anatomical differences between a group of patients and a control group. Several image preprocessing steps are included in the method's pipeline and deformable image registration is the one playing a crucial role. Its central idea is to find local forces which will deform a floating image to make it more similar to a reference image. The involved non-linear transforms are either based on smooth basis functions [3], [4], [5] or they are physically interpreted, e.g. by mechanics of continuum [6], [7], [8]. The former group of methods produces smooth low-dimensional deformations which are able to suppress only gross anatomical inter-subject differences, whereas the goal of the latter methods is to achieve a perfect match. In [9], [10], images in VBM are put into a stereotaxic space by an affine transform and then they are warped to the reference image by low-dimensional parametric deformations based on lowest-frequency components of discrete cosine transform. The coefficients are searched in an optimization algorithm which minimizes the residual squared difference between the images and simultaneously maximizes the smoothness of the deformations. Only one scaling parameter is incorporated to count for differences in intensities of the images, what makes it suitable for monomodal images only.

In this paper, we propose a deformable registration algorithm proper for multimodal images. Below, we first explain the methods used in our algorithm and then we present experimental results obtained from its evaluation.

Methods

Our deformable registration is performed by a multilevel block-matching technique, see fig 1. A floating image N is deformed to match a reference image M in an iterative process. A resulting displacement u is made up from local translations of the blocks of the floating image N by radial basis function (RBF) interpolation. The translations representing warping forces f are found by maximizing symmetric regional similarity measures. The floating image N is assumed to be brought into the same coordinate space as the reference image M by a previous linear registration step.

Regional symmetric similarity measures

where S_H is a point similarity measure derived from the global joint entropy of the images. All the defined point similarity measures depend on the joint intensity distribution, which is estimated from the joint histogram, which is not known until the images are perfectly matched. Thus, it is usually estimated from the images aligned by the previous registration step. In this way, the deformable registration is done also in [12], where a region similarity measure based on conditional probabilities is proposed. It is rewritten here as another point similarity measure:

At each level of subdivision, translations of rectangular blocks of the floating image N are searched in an optimization algorithm, which maximizes a selected region similarity measure. Inspired by symmetric registration proposed in [13], the symmetric regional similarity measure is obtained here as a sum of two partial similarity measures. These are computed in the blocks of the floating image according to the reference image blocks as well as in the reverse direction.

To avoid getting trapped in local minima, a combination of extensive search and hillclimbing algorithms is used here. First, a space of all possible translations is searched with a relatively large step. P best points are then used as start points for the following hillclimbing. The minimum of P local minima obtained by the hillclimbing is then declared as the global minimum.

Multilevel deformation

where u_k(x) is the displacement of a grid point x in the k^th dimension, R is the radial basis function of the distance ||x-w_i|| between the grid point x and the center of the i^th block w_i. The coefficients α_i are computed by putting the translations f into (6) and solving the resulting linear system of B equations separately for each dimension k. The compactly supported Wendland's RBF, which was successfully used for landmark-based deformable registration in [4] is used here. Its mathematical properties hold for different spatial support, which is important for the multilevel strategy. For each level of subdivison, the block size is set to the half of the size at the previous level. The displacements are gradually incremented over all levels, refining the resulting deformation in the coarse-to-fine manner. The regions containing poor contour or surface information can be eliminated from the matching process and the algorithm can be accelerated in this way. The subdivision is performed only if at least one voxel in the current region has its normalized gradient image intensity bigger then a certain threshold, see fig. 2.

Tissue probability maps

The computation of the global joint histogram is not the only way how to estimate the joint intensity distribution of the images M and N. When the registration is done in a stereotaxic space, tissue probability maps are often available, representing a kind of prior information, which can be used. Here, another estimate of the joint intensity distribution is made with their use. It is further combined with the usual estimate obtained from the global joint histogram by a weighting parameter λ:

The joint intensity estimates are re-computed in the each iteration of the registration algorithm. As the intensities of the floating image are not spread on a regular grid, simple unitary increments of individual histogram bins cannot be done. Instead of that, all histogram bins corresponding to the reference image voxels in the neighborhood of a displaced floating image voxel are increased by a value equal to the value of an interpolation kernel function. Cubic B-splines are used for that purpose in the generalized partial volume interpolation (GPVE) described in [14]. This approach is used here to compute the global joint histogram. In the computation of the joint distribution estimate based on the tissue probability maps, the interpolation has to be done among the voxels of the floating image. Thus, some ideas from the distance-weighted scattered data interpolation methods were adapted here and the locally bounded kernel function described in [15] is used in the partial volume interpolation scheme.

Experiment and results

The performance of the proposed algorithm with various similarity measures was evaluated with the use of 2D image data obtained from Simulated Brain Database (SBD) [16]. The original size of the SBD transversal slices is 181×217 pixels with the pixel size 1×1 mm. For the evaluation purposes, the images were padded to the size of 217×217 pixels. The square image size is convenient for the subdivision scheme. Synthetic deformations were composed from random translations at 10% randomly selected pixels in transversal slices. These force fields were smoothed by gaussian filtering with random standard deviation (σ=10±5 mm) to obtain final displacements which were then applied on 20 intensity images and corresponding segmented images. The average initial overlap error was 41.0%. The deformations were then recovered by the proposed deformable registration algorithm and the average decrease of the overlap error Δe was computed. T1 weighted, T2-weighted and PD (proton density) images with 3% noise and 20% intensity nonuniformity were used as floating images and the T1-weighted image with no noise and no intensity nonuniformity was used as the reference image. The results for various levels of decomposition are summarized in tab 1. The overlap error got smaller up to the 5^th level, when the subimages had size of 7x7 pixels. Although the next level gave an increase in the global mutual information, the alignment measured by overlap errors and also by visual inspection was constant or worse. This level was considered as the maximum subdivision level for this algorithm.

Tab. 1. Quantitative validation of the proposed algorithm. The averages of the overlap error decrease Δe were computed for various similarity measures, various image pairs and various levels of decomposition.

Conclusion

An algorithm for low-dimensional atlas-based registration of MRI images was presented. Four various symmetric region similarity measures were studied in an experiment in which synthetic deformations were recovered and the performance of the algorithm was quantified by the decrease of overlap error in segmented images. The similarities were measured with the use of joint histogram and tissue probability maps from MRI brain atlases. The overlap error was lower with the similarity measures S_PC and S_PMI depending on probabilities than when the similarity measures S_MI and S_UH depending on uncertainty were used. In our implementation, partial volume interpolation scheme was used, so that it was unnecessary to compute the deformed floating image during the registration process. The proposed algorithm is suitable for the voxel based morphometry, as the precision of the registration can be controlled by the maximum level of decomposition. Thus, only gross inter-subject anatomical differences can be suppressed and the important variability for statistical parametric tests can be preserved. In addition, the use of multimodal similarity measure allows to use an arbitrary available brain atlas without any need to transform intensities in the images to obtain monomodal data.

Acknowledgement

This work has been partly supported by the grant projects 1ET2085120511 AC CR, 102/04/0472 and 305/04/1385 from GACR, and the Research Programme of Brno University of Technology MSM 0021630513.

Low-dimensional Multimodal Deformable Registration of MRI Brain Images in Stereotaxic Space

Introduction