Global image registration using a symmetric block-matching approach

Abstract. Most medical image registration algorithms suffer from a directionality bias that has been shown to largely impact subsequent analyses. Several approaches have been proposed in the literature to address this bias in the context of nonlinear registration, but little work has been done for global registration. We propose a symmetric approach based on a block-matching technique and least-trimmed square regression. The proposed method is suitable for multimodal registration and is robust to outliers in the input images. The symmetric framework is compared with the original asymmetric block-matching technique and is shown to outperform it in terms of accuracy and robustness. The methodology presented in this article has been made available to the community as part of the NiftyReg open-source package.


Introduction
Medical image registration is core to many image analysis pipelines.It consists of bringing two or more images into spatial alignment, often mapping one image into the space of another.However, recent studies have highlighted that this directionality in the registration can create a bias in subsequent analyses.
Yushkevich et al. 1 and Leung et al. 2 highlighted the effect of asymmetric processing when estimating the brain atrophy occurring during a neurodegenerative process such as Alzheimer's disease.Thompson et al. 3 pointed out that bias in registration directionality has a large impact on the extraction of imaging biomarkers derived from tensor-based morphometry.Fox et al. 4 enumerated several assessment methods, including symmetry, against which new methodologies for atrophy estimation should be tested.
Several approaches have been proposed to address this issue in the context of nonlinear registration.Christensen and Johnson 5 jointly optimize forward and backward transformations while minimizing an inverse-consistency error criterion.Tao et al. 6 and Modat et al. 7 used similar schemes using the demons approach 8 and a cubic b-spline parametrization, 9 respectively.Vercauteren et al. 10 took advantage of a nonstationary velocity field parametrization to concurrently optimize forward and backward transformations through the exponentiation of a velocity field and its negated version.Avants et al. 11 presented a symmetric approach based on the optimization of a nonstationary velocity field discretized into intermediate space images.
Little work has, however, been proposed to remove directionality bias in the case of global registration.Leung et al. 2 proposed to run several asymmetric registrations, for example forward and backward in the case of pairwise registration, and to combine the obtained transformations into a consensus registration through averaging in the log-Euclidean space.Reuter et al. 12 developed a global registration technique based on robust statistic optimization performed in a midpoint space.This approach minimizes the residual differences between input images while rejecting outliers.A linear intensity scaling is used in order to increase the robustness of the algorithm and deal with different ranges of intensity within the same modality or pulse-sequence.This method has been altered by Wachinger et al. 13 to accomodate for multimodal images by using local entropy.
In this paper, we propose a robust and symmetric registration method that differs from Reuter et al. in two main aspects.First, we used a block-matching approach 14 to establish the spatial correspondences, where the normalized cross correlation is used as a measure of similarity.Due to the small dimension of the blocks under consideration, the proposed approach is suitable for multimodal registration applications. 15The second main difference is that joint forward and backward transformation parameters are simultaneously calculated rather than using a midpoint.This removes the need to discretize the transformed input images into an average space, which can be problematic when input images have different fields-of-view and resolution, as often seen in registration of different modalities.In these scenarios, the user must make an algorithmic decision on how to discretize the midspace, which can lead to a bias toward one of the images.Assuming two images of different resolutions, the user has to decide how to discretize the midspace.With the proposed approach, the transformation is optimized within the native spaces of the input images.Similarly to Reuter et al. and Wachinger et al., the proposed approach is also robust to the presence of outliers thanks to the use of a least-trimmed squared (LTS) optimization.
We evaluated the accuracy, precision, and robustness of our symmetric block-matching-based approach using several databases of synthetic magnetic resonance (MR) images, MR longitudinal studies, multimodal studies of magnetic resonance imaging (MRI), Positron emission tomography (PET), and computerised tomography (CT) scans, and pairs of preoperative and intraoperative MR.

Method
2.1 Block Matching for Global Registration: the Classical Approach In the original block-matching method, the registration is performed using a two-step approach.The first step is to establish point correspondences between the current transformed floating image and the reference image using the block-matching strategy.This is done by dividing the reference and transformed floating images into blocks of uniform size.Each block in the reference image is compared with image blocks in the corresponding neighborhood of the transformed floating image.The matching block is the one with the maximum normalized cross correlation (NCC) with the reference block computed as: where b r and b f denote blocks in the reference and warped floating image, respectively, μ and σ correspond to the mean and standard deviation value within a block, and N is the number of voxel in a block.The second step is to compute the transformation parameters from these point correspondences using a LTS regression method, which ensures robust outlier rejection.This optimization technique is also computationally efficient as each LTS iteration has an analytical solution.These two steps are iteratively performed, where the optimization starts with a large block search neighborhood (corresponding to gross displacements) and then moves toward a smaller search neighborhood (corresponding to finer displacements) in a coarse to fine strategy.

Directionality Bias in Image Registration
The classic block-matching image registration maps one image into the space of another, producing a transformation T I→J that maps image I to image J and enables image J to be warped into the space of image I.If the two images were reversed so that J is mapped to I, a new transformation T J→I is produced.The directionality that is prevalent in many image registration algorithms would result in T I→J ≠ T −1 J→I .If an algorithm is symmetric, then these two transformations would be the inverse of each other, and the registration result would not depend on the order of images.
Figure 1 illustrates the impact of the directionality bias on a subsequent analysis: the full brain atrophy estimation.In this example, two T1-weighted MRIs from the same subject and acquired with a 1-year interval are spatially normalized to a midspace based on an affine transformation.This step ensures that no resampling bias affects the brain atrophy estimation.Three registration approaches are used to estimate the affine parametrization: (i) an asymmetric registration where the baseline scan is considered as the reference and the follow-up scan as the floating, T asym;b→f , (ii) an asymmetric registration where the follow-up scan is considered as the reference and the baseline scan as the floating, T asym;f→b , and (iii) a symmetric registration scheme, T sym .After the resampling of the input images to a midspace, they are corrected for differential bias fields and the atrophy is estimated using the boundary shift integral method as described by Leung et al. 2 The obtained full brain volume changes over one year are 14.49, 20.55, and 16.51 mm using the forward, backward, and symmetric registration schemes, respectively.It can be noted that using the asymmetric registration approaches, the obtained atrophy measurements are different even though the composition of the obtained matrices for this specific case is close to identity, hence close to being inverse consistent 1.0004 0.0004 0.0006 0.0192 0.0004 0.9944 0.0011 −0.0516 0.0002 0.0016 0.9969 0.0437 0.0000 0.0000 0.0000 1.0000 where ∘ denotes the composition operator.

Symmetric Extension to Block-Matching Global Registration
Since the original block-matching-based registration introduces a directionality, the results would yield transformations that are not the inverse of each other when the order of images are switched.The inverse consistent extension to the block-matching algorithm ensures that the registration is symmetric and no bias is coming from the registration direction.Similar to the classical approach, the optimum transformation is estimated in an iterative fashion where the current estimated transformation is updated at each iteration.Each iteration follows a two-step scheme where the first step is to establish Fig. 1 Two scans from the same subject acquired with a 1-year interval are affinely registered using two asymmetric and a symmetric approaches.Using the obtained transformations, they are spatially normalized to midspaces before computing the full brain volume changes using the boundary shift integral method.
point correspondences between the two images.Using block matching, two sets of correspondences are estimated: f CI→J g mapping points from image I to image J and f CJ→I g to map points from image J to image I.This two sets of correspondences are assessed from blocks defined in image I and J, respectively.
The second step is to update the transformation parameters estimated through LTS regression.In the original block-matching-based approach, the updated parameters are incremental and composed with the current estimate of the transformation.For example, if the current transformation at iteration t is F t and an update f t is computed, then the updated transformation would be F tþ1 ¼ f t ∘ F t .If the same was applied to the inverse of the transformation B t ¼ F −1 t , and an incremental update b t was obtained such that b t ¼ f −1 t , then the properties of matrix multiplication state that This results in the incremental updates being premultiplied in one direction but postmultiplied in the other direction.Since matrix multiplication is not commutative, this approach would break symmetry as the resulting transformations would depend on how the images were ordered.To address this issue, the block-matching correspondences for the symmetric approach are always established using the original image positions at every iteration.This is achieved by combining, at every iteration, the obtained block-matching correspondences with the previous estimate of the transformation.As a result no matrix multiplication has to be performed since the update is integrated into the LTS regression as: instead of the previous update approach: To ensure inverse consistency between the transformation parameters, the transformation matrices obtained through the LTS fitting are averaged 16 and updated as where Expm and Logm are the exponential and logarithm matrix operators, respectively.Using such an approach, the result transformations, forward and backward, are inverse consistent (T I→J ¼ T −1 J→I ) and the directionality of the registration does not affect the recovered transformation parameters.

Validation
The proposed block-matching method has been implemented as part of the NiftyReg open-source software. 17All registrations were run using the same parameters to allow for a fair comparison of the results.NCC was used as a local measure of similarity in the block-matching procedure.The block dimensions were chosen to be 4 × 4 × 4 voxels.The blocks were sorted according to decreasing intensity variance and the top 50% of the blocks in the reference images were considered for the LTS optimisation step.During the LTS optimization, 50% of the blocks with the largest-squared Euclidean residuals (according to the current estimated transformation) were considered as outliers.For the registrations, a multiscale pyramidal approach with three resolution levels was used to avoid local minima and increase the capture range.The registrations produced a full affine transformation (12 degrees-of-freedom).For the retrospective image registration evaluation (RIRE) database, 4 resolution levels and a rigid transformation (6 degrees-of-freedom) were used.

Validation on Synthetic Data
In order to assess the capture range and robustness of the proposed symmetric approach, two synthetic images 18 were used: a T1-weighted (T1w) and a T2-weighted (T2w) MRI.Both images were simulated with 9% noise level and 20% intensity nonuniformity.Orthogonal views of these images are shown in Fig. 2. The images were a priori in alignment and known affine transformations were applied to one image.Registrations were then performed using the symmetric and the asymmetric approaches in both forward and backward directions.The registration error was computed at the eight corners and the average Euclidean distance between the ground truth and the recovered positions for these corner points are reported.
The known transformations were generated by applying rotations from −45 to 45 at 15 increment steps, translations from −45 to 45 voxels with a 15 voxels increment, scaling from 50% to 150% with a 10% increment, and finally shearing from −0.2 to 0.2 with a 0.05 increment.This resulted in 4851 transformations (varying each of the 12 components independently) and a total of 14,553 individual registrations (symmetric, asymmetric forward, and asymmetric backward).The two input images were also registered using the FMRIB's linear image registration tool (FLIRT) algorithm, 19 part of the FSL software package. 20The normalized mutual information was used as a similarity measure.Note that due to the high level of noise and intensity nonuniformity, FLIRT lead to an error of 1.30 and 1.16 voxels for the forward and backward approaches, respectively, in the case where the input images were already aligned.Based on these error values, any registrations results with an average error larger than 2 voxels were classified as a failure.The success rates for the symmetric, forward, and backward implementations were 91.12%, 78.96%, and 86.02%, respectively.For 1742 transformations, the symmetric approach successfully recovered the transformation when one asymmetric scheme failed.In 410 cases, the symmetric recovered the known transformation when both the forward and backward scheme failed.For 19 occurrences, the symmetric approach failed whereas the backward method successfully recovered the transformation and, in one case, the symmetric registration failed while both asymmetric approaches recovered the transformation.

Validation on the RIRE Database
The RIRE project, based on the Vanderbilt database, 21 consists of CT, PET, and MRI scans with an associated marker-based gold standard transformation. 22We assessed the accuracy of the proposed symmetric approach using this database and compared its result with the asymmetric approach.For all registrations, the CT and PET images were used as references and the landmark errors were assessed in their space.The voxel dimensions of the CT and PET images were 0.65 × 0.65 × 4.00 and 2.59 × 2.59 × 8.00 mm, respectively.The proton density, T1w, and T2w MRI were used as floating images.Figure 3 shows coronal views of multiple modality scans from a subject.
Intrasubject registration was performed on 9 subjects and the registration error, defined as the Euclidean distance in millimeters, between the ground truth and the recovered transformed positions, was computed at 704 landmark positions.Table 1 shows, for all registration techniques, the mean and maximum errors per modality and the overall error.Figure 4(a) shows the cumulative registration error over 364 landmarks when the CT images are used as reference, and Fig. 4(b) presents the same result for 340 landmarks defined on the PET images.

Validation on the MIRIAD Database
We used the MIRIAD database, 23,24 to assess the transitivity properties of the symmetric and asymmetric approaches.The MIRIAD database consist of brain MRI scans from 46 subjects diagnosed with Alzheimer's disease and 23 age-matched healthy controls.For all subjects, we used a baseline, 6-and 12-months followup scans.The 6-and 12-months follow-up scans were registered to their corresponding baseline (T BL−6 and T BL−12 respectively) and the 12-months scans were registered to the 6-months scans (T 6−12 ).Transitivity error was defined as the average Euclidean distance at the transformed corners of the baseline image using the direct correspondences (T BL−12 ) compared with the composed transformation (T 6−12 ∘ T BL−6 ).The mean transitivity error over the 69 subjects was 0.54 (0.27) and 0.66 mm (0.34) for the symmetric and asymmetric approaches.The symmetric approach leads to significantly lower transitivity error than the asymmetric approach (two-tailed, paired t-test, p < 0.005).The minimum and maximum average errors were [0.15, 1.35] and [0.20, 1.60] for the symmetric block matching and the asymmetric block matching, respectively.For comparison, using FLIRT with the default parameters, we obtained a mean error of 1.00 mm (0.71) and minimum and maximum average errors of [0.27, 4.74].

Intraoperative Images
A dataset comprising 15 pairs of pre and intraoperative T1weighted MR images was used to assess the robustness of the proposed algorithm to missing tissues.The preoperative MRI were acquired on a 3T GE Signa Excite HD (General Electric, Waukesha, Milwaukee, Wisconsin) with a spatial resolution of 0.9 × 0.9 × 1.1 mm.The intraoperative scans, acquired on a 1.5T Siemens Espree (Erlangen, Germany), have a spatial resolution of 1.1 × 1.1 × 1.3 mm 3 and were all acquired after some tissue resection was performed.Note that the intraoperative acquisitions do not have the brain center in the field-ofview.Figure 5 shows the pre and intraoperative scans from one subject.
We performed the 15 intrasubject registrations using the proposed symmetric block-matching approach and twice with the asymmetric block-matching approach (one where the  preoperative image served as reference and one where the inraoperative was the reference).It is difficult to localize anatomical landmarks in the intraoperative MRI images due to severe nonlinear geometric distortions and degraded image quality, hence we only performed a qualitative analysis of the registration results.Two reviewers were asked to qualify the registration results as successful, acceptable or failure based on visual assessment while being blinded to the registration scheme.
The results are reported in Table 2.The symmetric was able to successfully register all 15 datasets, while multiple cases using the asymmetric registrations were classified as failures.

Discussion
We presented a symmetric extension of the block-matching algorithm for global registration.To the best of our knowledge, the proposed scheme is one of the two available implementations for robust multimodal symmetric registration, the other being the implementation of the work by Wachinger et al. 13 The algorithm has been compared with the original asymmetric framework in both mono and multimodal registration.Using simulated transformations, we found that the symmetric formulation led to an increase in capture range.This method also yield a lower transitivity error over multiple time-points and multiple subjects, a small but acknowledged source of bias in subsequent analyses.Using the RIRE database for evaluation of intermodality registration, we found that the symmetric approach outperformed its asymmetric counterpart both in terms of accuracy and in term of robustness as shown by the lower cumulative errors   over multiple landmarks.Finally, using pre and intraoperative MRI data, the proposed algorithm was shown to be robust to missing tissue.Note that all validations were performed without using any masking or initial alignment of the brain center of mass, as the overall aim was to compare several methods in the same context.As a result, we do not claim that the proposed implementation is suitable for all applications and outperforms all other implementations.Our aim was to highlight the advantages of using a symmetric registration approach as opposed to an asymmetric approach.

Fig. 3
Fig. 3 Corresponding coronal slices of the same subject in the multiple modalities available as part of the retrospective image registration evaluation database.

Fig. 4
Fig. 4 Error in millimeters over multiple landmarks using the symmetric and asymmetric block-matching approaches.(a) and (b) The errors assessed when the images used as reference are CT and the PET, respectively.

Table 1
Registration errors from multimodal intrasubject registrations using the RIRE database.Sym.BM: symmetric block-matching registration, Asym.BM: original, asymmetric block matching.The columns show the mean and maximum errors of the various multimodal combinations registered, with the overall error in the last column.

Table 2
Qualitative evaluation of pre and intraoperative registration performance.