Dual-core steered non-rigid registration for multi-modal images via bi-directional image synthesis

1.1.1. Information theory based registration methods

The first category is the information theory based methods, which use mutual information (MI) as the similarity metric for multi-modal image registration (Collignon et al., 1995; Wells et al., 1996; Maes et al., 1997; Viola and Wells III, 1997). Based on this similarity metric, rigid and affine registration results can be efficiently obtained between multi-modal images. However, MI is a good global similarity measurement, which has limited power in capturing local anatomical details, and thus has limited capacity in tackling local deformations. To address this issue, normalized mutual information (NMI) (Studholme et al., 1999) and localized mutual information (LMI) (Klein et al., 2008) have been proposed to measure the similarity between local anatomical structures. However, insufficient number of voxels within the image patch undermines the estimation of intensity distribution when computing MI. As a result, local anatomical structures are often not well aligned by MI-based registration methods. A good review of MI-based image registration methods can be found in (Pluim et al., 2003).

Besides MI, Chung et al. (2002) proposed using Kullback–Leibler distance (KLD) as the similarity metric to guide multi-modal image registration, by enforcing the joint intensity distribution of the to-be-aligned image pair following the priors of learned distributions from the well-aligned images. In addition, So and Chung (2011) further proposed using Bhattacharyya distance (BD) as the similarity metric, which obtained better registration performance compared with KLD-based registration.

However, either MI-based or KLD-based (BD-based) registration methods have a common limitation that these methods are all based on the intensity probability distribution, and ignore spatial information of anatomical structures. This means that different anatomical regions will be found as being matched as long as they have similar intensity distributions.

1.1.2. Image synthesis based registration method

For this kind of methods, one image modality is synthesized from the other modality in order to reduce the large appearance gap between different modalities. Afterwards, the multi-modal image registration problem is simplified to a mono-modal registration problem, where most existing registration methods can then be applied. Roche et al. (2001) synthesized an ultrasound image and then estimated rigid registration between the ultrasound and the MR images. Wein et al. (2008) simulated the ultrasound image from the CT image to solve the rigid and affine registration problems.

There are two main problems with the current image synthesis based registration methods. The first problem is that the image synthesis is always performed in a single direction. This will often introduce bias when performing registration only based on one modality image, by ignoring all the anatomical details in the other modality. The second problem is that most of the existing works focus on rigid/affine registration. Few work tackles the non-rigid registration problem. Since existing methods often synthesize from the image modality with rich anatomical details (e.g., MRI) to the image modality with limited anatomical details (e.g., ultrasound/CT), it is quite difficult to accurately estimate local deformations based on the image modality with limited anatomical details, especially for the case of pelvic images where the main pelvic organs often have large local deformations. Therefore, an accurate and robust image synthesis method, which can preserve adequate anatomical details, is essential for image synthesis based multi-modal non-rigid registration.

1.1.3. Image synthesis methods

To date, many image synthesis methods are proposed although most of them are not for registration purpose.

1.1.4. Other multi-modal registration methods

Some feature-based registration methods (Shen and Davatzikos, 2002; Zacharaki et al., 2009; Ou et al., 2011; Heinrich et al., 2012; Toews et al., 2013) estimate deformation fields by leveraging different artificial features (e.g., geometric moments (Shen and Davatzikos, 2002), Gabor attributes (Ou et al., 2011)), rather than directly using image intensities. In this way, the influence of appearance differences between modalities can be minimized during registration. However, it is difficult to find effective features that can describe common anatomical details and are also invariant to the appearance differences. Additionally, some methods conduct the registration and segmentation jointly (Ashburner and Friston, 1997; Jia et al., 2012), or use the segmentation images (Shen and Davatzikos, 2002) to help improve the multi-modal registration. But, accurate segmentation is needed for accurate non-rigid registration, which is often not easy to achieve.

2.2.1. Patch-wise random forest (P-RF) regression

The basic theory of random forest can be found in (Liaw and Wiener, 2002). It can be used for non-linear classification and regression, and has been widely adopted in disease diagnose (Zhang et al., 2016), image registration (Wei et al., 2016) and segmentation (Wang et al., 2015; Gao et al., 2016). Here, we directly introduce the improved patch-wise random forest for the application of image synthesis. Overall, in MR image synthesis from CT, it can be used to regress the MRI patch intensity from the corresponding CT patch and build a nonlinear mapping to bridge the appearance gap between the two modalities based on multiple decision trees. Fig. 3 illustrates the training and testing stages of patch-wise random forest for the MR image synthesis from CT.

In the training step of patch-wise random forest, the CT and MR images of the same patient are already pre-aligned (with the details of the CT and MRI pre-alignment are described in Section 3.2). Based on this pre-aligned dataset, we first randomly select N voxels to obtain N training samples. Each sample consists of two components: 1) the input appearance feature vector extracted from a single CT patch, which can be denoted by x, and 2) the target MRI patch intensity value corresponding to the center of the CT patch, which can be represented by vector y. For training the non-linear mapping from CT to MRI, the input is N CT feature vectors X = [x₁, x₂, ⋯, x_N], along with the corresponding N target MRI patch intensity values Y = [y₁, y₂, ⋯, y_N]. Random forest consists of multiple decision trees, with each one trained independently. For each tree, its training is conducted by learning a set of splitting nodes to recursively partition the training sample set. Specifically, in each split node, for a feature indexed by k, its optimal threshold τ is found to best split the training set into left and right subsets S_L and S_R with consistent target MRI patch intensity values. Mathematically, it is to maximize the variance reduction by a split:

where V(·) computes the averaged variance of target MRI intensity values across each element of the patch vector in the training set, x^k indicates the k-th feature, and S indicates a training set. Once the optimal combination of feature index and its threshold is found for a split node, the training set S is partitioned into S_L and S_R with the sample numbers N_L and N_R, respectively. The same split operation is recursively conducted on S_L and S_R, until 1) the tree reaches the maximum tree depth, or 2) the number of remaining training samples is too few to be split.

In the testing step, given a testing sample with feature vector x_new extracted from a new CT patch sample (see the right part of Fig. 3), it is pushed to the root split node of each tree in the forest. Under the guidance of the split node (i.e., go left if ${\mathit{x}}_{\text{new}}^k<\tau$, and go right otherwise), the testing sample will arrive at a leaf node of each tree, where the averaged target MRI patch intensity values of training samples in that leaf are used as a prediction of the tree. The final output of random forest is the average of patch predictions from all trees. Obviously, by going through all the voxels in the new CT image, the predicted patches are highly overlapping. Thus, by averaging all those overlapping values, the predicted MRI can be finally obtained.

Since the target MR image has strong spatial clues, we combine the patch-wise representation with random forest for the MR image prediction. The advantage of patch-wise random forest regression is that, by predicting a MRI patch as a whole, the neighborhood information can be naturally preserved and lead to better image synthesis eventually. Moreover, since more target values can be obtained as the output for each tree, the number of trees in the whole forest can be reduced, which improves the efficiency of the prediction.

2.2.2. Feature extraction

Feature extraction is a crucial step for random forest training. An effective feature should have the capacity to distinguish different anatomical structures, as well as robustness to noise. Typically, it is not effective to directly use patch intensities as appearance features, since soft tissues have low contrast in the CT image, and also the CT image is always noisy. Therefore, the individual CT values are insufficient to characterize different pelvic structures, as they are also sensitive to noise. To solve this problem, we extract Haar-like features from CT patch to serve as appearance features for random forest training. Specifically, a Haar-like feature (Gao et al., 2016) describes 1) the averaged intensity within a sub-block in the patch, or 2) the averaged intensity difference between two sub-blocks in the patch, as shown in Fig. 4.

By using the integral image (Viola and Jones, 2004), in which each position contains the summation of the voxel intensity values from the original position to the current location, Haar-like feature can be calculated with the constant time complexity. To generate more Haar-like features, we randomly sample the information within the patch. Since random forest has the inherent capability of feature selection, only the meaningful Haar-like features will be retained. It is worth noting that, to cover both local and global appearances for the underlying voxel, Haar-like features are extracted from coarse, medium and fine image resolutions by using the same patch size, respectively.

2.2.3. Auto-context model (ACM) based image synthesis refinement

In the patch-wise random forest, the MR image patch of each voxel is estimated from the Haar-like features extracted from the corresponding CT image patch. To build better non-linear mapping between the CT patch features and the MR image patch, we further incorporate the neighboring prediction results to enhance the patch-wise random forest training. In particular, an auto-context model (Tu and Bai, 2010) is adopted to refine the synthesized MR image iteratively. In this paper, we perform three layers as illustrated in Fig. 5. Note that, the auto-context model is used only in the last two layers, as indicated by a dashed box.

In the first layer, the appearance features (i.e., Haar-like features) from CT images are extracted to train a patch-wise random forest. Then, the trained forest is used to generate the initial synthesized MR images. In the second layer, to coordinate the prediction information in the neighborhood, the Haar-like features (namely the context features, as shown in Fig. 5) are further extracted from those initial synthesized MR images, similar to the extraction of appearance features. By combining these context features with the CT appearance features, a second patch-wise random forest can be trained, which also leads to the update of the synthesized MR images and their context features. This process iterates until it reaches the maximum number of layers.

2.4.1. Intensity-based non-rigid registration

After the CT and MR images are synthesized, we have a total of four images: actual CT, actual MRI, synthesized CT from MRI, and the synthesized MRI from CT. Based on the actual and synthesized images, we can utilize the existing non-rigid registration methods to estimate the deformations 1) from the synthesized CT to the actual CT, and 2) from the actual MRI to the synthesized MRI. Among various existing registration methods, we choose two popular non-rigid registration methods for evaluation: 1) Diffeomorphic Demons (D. Demons) (Vercauteren et al., 2009), which is a non-parametric registration method widely used in various applications; 2) Symmetric Normalization (SyN) (Avants et al., 2008), which is a top-ranked method among the 14 state-of-the-art non-rigid registration methods evaluated for brain image registration in (Klein et al., 2009). For SyN, we use the Advanced Neuroimaging Tool (ANTs) (Avants, Tustison et al., 2011) for its implementation.

2.4.2. Dual-core deformation fusion (DDF) for MRI-to-CT registration

To overcome the limitation of multi-modal image registration methods that utilize only single-directional image synthesis, we utilize both synthesized CT and synthesized MRI in MRI-to-CT image registration, which can take the advantage of complementary information from both modalities. Let I_CT, Î_CT, I_MR and Î_MR denote the actual CT, synthesized CT, actual MRI, and synthesized MRI, respectively. Traditionally, a general objective function for the conventional MRI-to-CT image registration can be given as:

where φ is the deformation field to be estimated, ℳ is the dissimilarity metric, and φ(I_MR) means that deforms the subject MR image to the CT image space with the estimated deformation field φ. ℛ is a regularization term, used to constrain the smoothness of the estimated deformation field φ.

By using the bi-directional synthesized images, the deformation field can be estimated by simultaneously minimizing 1) the dissimilarity between the CT modality core and 2) the dissimilarity between the MRI modality core. Then, the registration objective function can be rewritten as:

To solve Eq. (4) and reuse the available registration tools, we apply the alternative optimization method. Specifically, we optimize the registration between the CT modality core (I _CT and Î_CT) and the registration between the MRI modality core (Î_MR and I_MR) separately as below:

Eqs. (5) and (6) are used to minimize the differences 1) between the CT image I_CT and the warped synthesized CT image φ_CT (Î_CT), and 2) between the synthesized MR image Î_MR and the warped MR image φ_MR (I_MR), respectively. In order to use the complementary information from both modalities, the dual-core deformation fusion is proposed accordingly:

Eq. (7) is used to ensure that the final deformation pathway φ to be close to both separately estimated φ_CT and φ_MR.

Both φ_CT and φ_MR can be optimized by using either D. Demons or ANTs-SyN, although the objective functions are slightly different in Eqs. (5) and (6). After estimating φ_CT and φ_MR, the final deformation φ can be effectively optimized by letting the gradient of Eq. (7) equal to zero, which brings to:

However, simple average of deformation pathways may make the well-aligned structures in one modality become misaligned after averaging. To alleviate this issue, we approximate the optimal solution of Eq. (4) by alternating the dual-core deformation fusion procedure (Eqs. (5)–(8)) based on the warped synthesized CT and the original MR images until convergence, as summarized in Algorithm 2.

For iteration i, the tentatively warped images ${\widehat{I}}_{\text{CT}}^{i-1}$ and $I_{\text{MR}}^{i-1}$ are used to estimate a next set of deformations ${\phi }_{\text{CT}}^i$ and ${\phi }_{\text{MR}}^i$. The estimated deformations are then merged to form a combined deformation φⁱ_, which is used to update the currently estimated deformation pathway φ = φ ∘ φⁱ. Here, “∘” denotes deformation field composition (Vercauteren et al., 2009), which means updating deformation field φ by φⁱ. Using this deformation field composition can help preserve the quality of the subject image, by warping the original moving image Î_CT and I_MR for just one time. The whole registration procedure iterates until the incremental deformation φⁱ is small enough.

3.2.1. Parameter setting

In image synthesis, the aligned CT and MR image pair is used to train our image synthesis models: 1) synthesizing CT from MRI (MRI-to-CT) and 2) synthesizing MRI from CT (CT-to-MRI). Table 2 provides parameters of patch-wise random forest (P-RF) training for bi-directional image synthesis. Specifically, 2-layer auto-context model (ACM) and 10-fold cross validation are applied. For P-RF training, the input feature sample is extracted from the local 15 × 15 × 15 patch and the output is the 3 × 3 × 3 patch intensity. Here, more trees and deeper tree depth are used for synthesizing MRI from CT, since this simple-to-complex mapping is more difficult. The voxel-wise random forest (V-RF) is also applied as the baseline method. Note that V-RF utilizes the same parameter setting as P-RF, but its output is a single voxel value, not the 3 × 3 × 3 patch.

3.2.2. Evaluation

Two measurements are used to quantitatively evaluate the image synthesis performance: 1) mean absolute error (MAE) and 2) peak signal-to-noise ratio (PSNR), which can be defined as:

where I and Î denote the original image and the synthesized image, and the image size represents by l, m, and n in three dimensions. Note that a high-quality synthesized image should have lower MAE and higher PSNR.

3.2.3. Bi-directional image synthesis results

Table 3 provides the quantitative performance of CT image synthesis from MRI, and Table 4 presents the quantitative performance of MR image synthesis from CT. From these two tables, we can observe that the patch-wise random forest is effective in building non-linear regression models to fill up the appearance gap between CT and MRI modalities. Furthermore, the image synthesis performance is improved by using more layers with the auto-context model. For both CT and MRI synthesis, the improvement achieved by adding the first-layer auto-context model is more significant than adding the second layer. This can qualitatively demonstrate the fast convergence of the auto-context model. In practice, two-layer auto-context model always leads to the convergence.

Figs. 6 and 7 provide detailed visualization of the synthesized CTs and the synthesized MRIs for one subject, along with the error maps between the actual and the synthesized images. The results are consistent with the performance in Tables 3 and 4. From Figs. 6 and 7, we can observe that the boundaries of both bone regions and soft tissue regions are clearer and less noisy by using P-RF, compared with results by the V-RF based image synthesis method. Moreover, the quality of the synthesized image is further improved by auto-context model refinement. The final synthesized image is quite similar to the actual image, by also preserving the smoothness and continuity. Note that, to achieve this promising performance, a relatively small number of trees (as shown in Table 2) is used to train our image synthesis model. This well demonstrates the effectiveness of the P-RF method, which allows obtaining adequate output values by using only a small number of trees.

It is reasonable that the performance of the synthesized CT image is much better than the synthesized MR image. This is because 1) the texture patterns are much more complex in the MR image than in the CT image, thus difficult to estimate, and 2) also the appearance in the MR image is more diverse than in the CT image. As shown in Figs. 6 and 7, the error maps of synthesized MRI are fuzzy, and also the synthesized MR images are smoother than the actual ones. However, the boundaries of pelvic tissues in the synthesized MR images are clear, and the appearance relationship is consistent between the synthesized MR image and the actual MR image, which provide crucial reference information for guiding the subsequent image registration.

3.3.1. Evaluation

After the MR image is registered to the CT image, we use two kinds of measurements to evaluate the registration performance: 1) Volume overlap. Specifically, Dice Similarity Coefficient (DSC) between manual labels on CT (L_CT) and aligned MRI (L_MR) is used as the main measurement to evaluate the registration performance. Besides, sensitivity (SEN) and positive predictive value (PPV) are also used to compare the proposed method with the existing methods. 2) Surface distance. Specifically, the distance measurement is used to account for the boundary discrepancies between the CT labels and the aligned MRI labels. Here, symmetric average surface distance (SASD) and Hausdorff distance (HAUS) are calculated to further demonstrate the registration performance. Eqs. (11)–(14) below provide the definition of all these measurements.

Here, L means a binary label mask from prostate, bladder and rectum. d(u, v) is the Euclidean distance between voxels u and v. The unit of the distance measurement is presented by millimeter (mm). A good registration result should have higher volume overlap and lower surface distance.

3.3.2. MRI-to-CT registration results

As described in Section 2.4, two state-of-the-art non-rigid registration methods (D. Demons and ANTs-SyN) are used in our experiment, with their detailed parameter settings shown in Table 5. For the preprocessed data, the CT and MRI of the same patient are already linearly aligned. The DSC values of prostate, bladder, and rectum after linear registration are 85.66% ± 5.65%, 89.29% ± 2.18% and 78.44% ± 3.54%, respectively, and we regard these results as the baseline.

Fig. 8 illustrates the MRI-to-CT registration results under different layers of ACM for image synthesis, and different dual-core deformation fusion (DDF) iterations for image registration, respectively. The mean DSC values of prostate, bladder, and rectum are evaluated. Here, the registration results in Fig. 8(a) use 3 iterations (3-iter) DDF, while the registration results in Fig. 8(b) use 2 layers (2-layer) ACM.

As shown in Fig. 8(a), more layers of ACM lead to better registration accuracy due to better quality of the synthesized images. Fig. 8(b) demonstrates that the DDF framework improves registration performance of both D. Demons and ANTs-SyN iteratively. In practice, we found that it is enough to have 2-layer ACM in image synthesis and 3 iterations (3-iter) in DDF.

Table 6 provides the mean and standard deviation of DSC of the total 20 subjects for the three crucial pelvic organs. It can be observed that, for D. Demons, which is originally not applicable for multi-modal image registration, can now work well by introducing the synthesized image. For ANTs-SyN, using MI as similarity metric obtains reasonable registration results on the original CT and MRI. However, as shown in Table 6, better performance can be obtained by using the synthesized image. This demonstrates that using synthesized image can enhance the performance of multi-modal image registration method. Moreover, the best performance is achieved by using our proposed dual-core deformation fusion algorithm. The consistently higher DSC and the lower standard deviation by our proposed method demonstrate both its robustness and accuracy in solving multi-modal image registration problem.

We further provide the mean SEN, PPV, SASD, and HUAS values of three pelvic organs (prostate, bladder and rectum) in Tables 7 and 8 to compare our proposed bi-directional image synthesis based registration method with single-directional image synthesis based registration method.

From Tables 7 and 8, we observe that our proposed bidirectional image synthesis based method consistently obtains higher SEN and PPV and also lower SASD and HAUS, compared with the single-directional image synthesis based registration methods. This demonstrates that, by using the complementary information from both modalities, the performance is significantly improved for multi-modal image registration. It is worth noting that, for all cases, the PPVs always have higher values than SENs. This phenomenon occurs when the organ contour is under-segmented compared with the ground truth. It is reasonable in our MRI-to-CT registration case, since the manual labels on MRI are always smaller than those on CT. Due to the low tissue contrast in CT, the manual labels on CT image always include more regions due to the difficulty of identifying the organ boundaries.

Fig. 9 provides the visualized registration results after ANTs-SyN registration for 3 subjects in the dataset. For each subject, the first row shows the actual CT and MR images and the corresponding synthesized MR and CT images, which are all used in our proposed iterative DDF registration framework. The second row provides the registration results, i.e., the warped MR image, by using the actual CT and MR image (CT & MRI), single-directional image synthesis (CT & S-CT, S-MRI & MRI), and our bi-directional image synthesis based registration method (Proposed). The visualized registration results are corresponding with those shown in Table 6. The yellow contours indicate CT’s ground-truth labels of the prostate, bladder and rectum. The red contours denote the labels of warped MR image after registration. From the overlaps between the ground-truth labels and the warped (estimated) labels, we can also observe that using single-directional synthesized image can partially improve the registration result, compared to the case of only using actual CT and MR images. Furthermore, it is obvious that our proposed method can better preserve the structures of pelvic organs during registration, and also achieves much higher registration performance.