Quantitative Error Prediction of Medical Image Registration using Regression Forests

TL;DR

This paper introduces a regression forest-based method for automatically predicting local registration errors in chest CT scans, aiding in quality control and evaluation of medical image registration techniques.

Contribution

It presents a novel automatic approach using regression forests to quantitatively predict registration errors in medical images, specifically applied to chest CT scans.

Findings

Mean absolute registration error: 1.07 mm (SPREAD), 1.76 mm (inter-database)

Classification accuracy: 90.7% (SPREAD), 75.4% (inter-database)

Enables automatic quality control in large-scale image analysis

Abstract

Predicting registration error can be useful for evaluation of registration procedures, which is important for the adoption of registration techniques in the clinic. In addition, quantitative error prediction can be helpful in improving the registration quality. The task of predicting registration error is demanding due to the lack of a ground truth in medical images. This paper proposes a new automatic method to predict the registration error in a quantitative manner, and is applied to chest CT scans. A random regression forest is utilized to predict the registration error locally. The forest is built with features related to the transformation model and features related to the dissimilarity after registration. The forest is trained and tested using manually annotated corresponding points between pairs of chest CT scans in two experiments: SPREAD (trained and tested on SPREAD) and inter-database (including three databases SPREAD, DIR-Lab-4DCT and DIR-Lab-COPDgene). The results show that the mean absolute errors of regression are 1.07 $\pm$ 1.86 and 1.76 $\pm$ 2.59 mm for the SPREAD and inter-database experiment, respectively. The overall accuracy of classification in three classes (correct, poor and wrong registration) is 90.7% and 75.4%, for SPREAD and inter-database respectively. The good performance of the proposed method enables important applications such as automatic quality control in large-scale image analysis.

Tables

Table 1. Table 1: An overview of the proposed features. Averages and maxima are taken over boxes of diameter [2, 5, 10, 15, 20, 25, 30, 35, 40] mm for the features: MIND MIND \mathrm{MIND} , std ⁡ 𝑻 std 𝑻 \operatorname{\mathrm{std}}\bm{T} , std ⁡ 𝑻 L std superscript 𝑻 L \operatorname{\mathrm{std}}\bm{T}^{\mathrm{L}} , CVH, ℰ ​ ( 𝑻 ) ℰ 𝑻 \mathcal{E}(\bm{T}) , ℰ ​ ( 𝑻 L ) ℰ superscript 𝑻 L \mathcal{E}(\bm{T}^{\mathrm{L}}) and Jac. Mutual information measures are calculated in boxes of [5, 10, 15, 20, 25, 30, 35, 40] mm. SID and GID are computed using Gaussian derivatives with standard deviations in the range [0.5, 1, 2, 4, 8, 16] mm.

Feature	$N_f$
$\mathrm{MIND}$	18	9 average boxes + 9 maxima boxes
MI	32	NMI, NMIS, PMI, PMIS calculated over 8 boxes
$\mathrm{std}𝑻$	18	9 average boxes + 9 maxima boxes
$\mathrm{std}{𝑻}^L$	18	9 average boxes + 9 maxima boxes
CVH	18	9 average boxes + 9 maxima boxes
$ℰ\left(𝑻\right)$	18	9 average boxes + 9 maxima boxes
$ℰ\left({𝑻}^L\right)$	18	9 average boxes + 9 maxima boxes
Jac	18	9 average boxes + 9 maxima boxes
NC	8	calculated over 8 boxes
SID&GID	12	calculated over 6 sigma’s

Table 2. Table 2: Distribution of the reference registration errors in each database, used during testing.

Database-iters	correct		poor		wrong		total
SPREAD 20	$848789$	$\left(84.1\%\right)$	$102837$	$\left(10.2\%\right)$	$58059$	$\left(5.8\%\right)$	$1009685$
SPREAD 100	$904796$	$\left(89.6\%\right)$	$66467$	$\left(6.6\%\right)$	$38422$	$\left(3.8\%\right)$	$1009685$
SPREAD 500	$925840$	$\left(91.7\%\right)$	$51910$	$\left(5.1\%\right)$	$31935$	$\left(3.2\%\right)$	$1009685$
SPREAD 2000	$935676$	$\left(92.7\%\right)$	$46170$	$\left(4.6\%\right)$	$27839$	$\left(2.8\%\right)$	$1009685$
SPREAD together	$3615101$	$\left(89.5\%\right)$	$267384$	$\left(6.6\%\right)$	$156255$	$\left(3.9\%\right)$	$4038740$
DIR-Lab-4DCT 20	$521481$	$\left(84.5\%\right)$	$71282$	$\left(11.5\%\right)$	$24543$	$\left(4.0\%\right)$	$617306$
DIR-Lab-4DCT 100	$540989$	$\left(87.6\%\right)$	$61131$	$\left(9.9\%\right)$	$15186$	$\left(2.5\%\right)$	$617306$
DIR-Lab-4DCT 500	$553757$	$\left(89.7\%\right)$	$53067$	$\left(8.6\%\right)$	$10482$	$\left(1.7\%\right)$	$617306$
DIR-Lab-4DCT 2000	$561909$	$\left(91.0\%\right)$	$46679$	$\left(7.6\%\right)$	$8718$	$\left(1.4\%\right)$	$617306$
DIR-Lab-4DCT together	$2178136$	$\left(88.2\%\right)$	$232159$	$\left(9.4\%\right)$	$58929$	$\left(2.4\%\right)$	$2469224$
DIR-Lab-COPD ANTsBSplineSyN	2643	$\left(88.1\%\right)$	184	$\left(6.1\%\right)$	173	$\left(5.8\%\right)$	3000
DIR-Lab-COPD elastix-advanced	2420	$\left(80.7\%\right)$	259	$\left(8.6\%\right)$	321	$\left(10.7\%\right)$	3000

Table 3. Table 3: Distribution of the reference registration errors, used during training.

Database	correct		poor		wrong		total
SPREAD together	$589854$	$\left(58.0\%\right)$	$270523$	$\left(26.6\%\right)$	$156881$	$\left(15.4\%\right)$	$1017258$
DIR-Lab-4DCT together	$328055$	$\left(53.0\%\right)$	$232499$	$\left(37.5\%\right)$	$58929$	$\left(9.5\%\right)$	$619483$

Table 4. Table 4: Regression results for single features on the SPREAD database. The columns indicate the number of features ( N f subscript 𝑁 𝑓 N_{f} ), the mean absolute error (MAE), the accuracy (Acc) and the F1 score. The sub-indices c, p and w correspond to correct [0,3), poor [3,6) and wrong [6, ∞ \infty ) mm classes, respectively.

	$N_f$	MAE	MAEc	MAEp	MAEw	Acc	F1c	F1p	F1w
$\mathrm{MIND}$	18	$1.10\pm 1.97$	$0.76\pm 0.72$	$1.59\pm 1.39$	$6.50\pm 5.88$	$89.8$	$94.9$	$34.1$	$83.0$
MI	32	$1.20\pm 1.88$	$0.89\pm 0.71$	$1.53\pm 1.14$	$6.30\pm 5.58$	$87.9$	$93.9$	$30.1$	$79.9$
$\mathrm{std}𝑻$	18	$1.59\pm 2.79$	$1.15\pm 1.78$	$2.98\pm 4.06$	$7.60\pm 6.12$	$85.5$	$92.7$	$22.4$	$64.4$
$\mathrm{std}{𝑻}^L$	18	$1.51\pm 2.40$	$1.11\pm 1.34$	$2.49\pm 3.05$	$7.32\pm 5.79$	$86.7$	$93.4$	$18.3$	$70.7$
CVH	18	$1.93\pm 3.29$	$1.49\pm 2.22$	$1.82\pm 2.00$	$9.80\pm 7.19$	$75.2$	$87.2$	$16.9$	$37.0$
$ℰ\left(𝑻\right)$	18	$2.00\pm 2.80$	$1.61\pm 1.76$	$2.18\pm 3.12$	$8.52\pm 6.48$	$69.8$	$82.8$	$17.0$	$43.5$
$ℰ\left({𝑻}^L\right)$	18	$1.68\pm 2.85$	$1.19\pm 1.71$	$3.19\pm 3.28$	$8.34\pm 6.74$	$84.4$	$92.6$	$11.7$	$54.8$
Jac	18	$2.15\pm 3.15$	$1.72\pm 1.90$	$1.91\pm 2.27$	$10.03\pm 6.97$	$68.2$	$83.7$	$13.0$	$31.4$
NC	8	$1.38\pm 2.89$	$0.90\pm 0.71$	$1.70\pm 1.68$	$9.41\pm 9.15$	$88.2$	$94.3$	$28.5$	$77.0$
SID&GID	12	$1.30\pm 2.12$	$0.94\pm 0.90$	$1.82\pm 1.63$	$6.95\pm 6.02$	$89.9$	$95.1$	$24.9$	$74.3$

Table 5. Table 5: Regression results for groups of features on the SPREAD database. The columns indicate the number of features ( N f subscript 𝑁 𝑓 N_{f} ), the mean absolute error (MAE), the accuracy (Acc) and the F1 score. The sub-indices c, p and w correspond to correct [0,3), poor [3,6) and wrong [6, ∞ \infty ) mm classes, respectively. MD, NN and LR stands for modality dependent, neural networks and linear regression, respectively.

	$N_f$	MAE	MAEc	MAEp	MAEw	Acc	F1c	F1p	F1w
Intensity	50	$1.09\pm 1.88$	$0.77\pm 0.68$	$1.49\pm 1.26$	$6.20\pm 5.68$	$90.3$	$95.1$	$35.7$	$83.6$
Registration	108	$1.32\pm 2.35$	$0.90\pm 1.04$	$2.10\pm 2.71$	$7.76\pm 6.01$	$90.0$	$95.1$	$31.5$	$78.4$
Combined	158	$1.07\pm 1.86$	$0.76\pm 0.65$	$1.47\pm 1.22$	$6.12\pm 5.64$	$90.7$	$95.4$	$38.1$	$84.4$
Combined-no pooling	8	$1.24\pm 2.22$	$0.85\pm 0.73$	$1.72\pm 1.64$	$7.39\pm 6.62$	$89.4$	$94.8$	$32.6$	$79.1$
Combined+MD	178	$1.07\pm 1.83$	$0.76\pm 0.65$	$1.46\pm 1.20$	$5.95\pm 5.59$	$90.7$	$95.4$	$38.3$	$84.5$
Combined (LR)	158	$1.86\pm 2.03$	$1.58\pm 1.34$	$2.47\pm 2.21$	$6.12\pm 4.97$	$77.3$	$87.3$	$17.0$	$67.6$
Combined (NN)	158	$1.13\pm 2.07$	$0.74\pm 0.70$	$1.81\pm 1.67$	$7.08\pm 5.88$	$89.8$	$95.0$	$31.2$	$79.6$

Table 6. Table 6: Leave one feature out results of SPREAD data. The columns indicate the number of features ( N f subscript 𝑁 𝑓 N_{f} ), the mean absolute error (MAE), the accuracy (Acc) and the F1 score. The sub-indices c, p and w correspond to correct [0,3), poor [3,6) and wrong [6, ∞ \infty ) mm classes, respectively.

	$N_f$	MAE	MAEc	MAEp	MAEw	Acc	F1c	F1p	F1w
Combined	158	$1.07\pm 1.86$	$0.76\pm 0.65$	$1.47\pm 1.22$	$6.12\pm 5.64$	$90.7$	$95.4$	$38.1$	$84.4$
$-\mathrm{MIND}$	140	$1.18\pm 1.96$	$0.83\pm 0.66$	$1.56\pm 1.50$	$6.70\pm 5.69$	$90.2$	$95.1$	$36.2$	$83.0$
$-$MI	126	$1.10\pm 1.98$	$0.75\pm 0.67$	$1.54\pm 1.30$	$6.66\pm 5.84$	$90.6$	$95.3$	$37.0$	$84.2$
$-\mathrm{std}𝑻$	140	$1.08\pm 1.86$	$0.76\pm 0.65$	$1.46\pm 1.18$	$6.14\pm 5.65$	$90.7$	$95.3$	$38.1$	$84.3$
$-\mathrm{std}{𝑻}^L$	140	$1.08\pm 1.89$	$0.76\pm 0.65$	$1.46\pm 1.22$	$6.21\pm 5.73$	$90.6$	$95.3$	$38.3$	$83.7$
$-$CVH	140	$1.07\pm 1.81$	$0.75\pm 0.65$	$1.46\pm 1.21$	$6.06\pm 5.98$	$90.7$	$95.4$	$38.4$	$84.3$
$-ℰ\left(𝑻\right)$	140	$1.07\pm 1.86$	$0.76\pm 0.65$	$1.46\pm 1.21$	$6.13\pm 5.64$	$90.7$	$95.4$	$38.2$	$84.5$
$-ℰ\left({𝑻}^L\right)$	140	$1.08\pm 1.85$	$0.76\pm 0.65$	$1.47\pm 1.22$	$6.12\pm 5.61$	$90.6$	$95.3$	$37.5$	$84.3$
$-$Jac	140	$1.08\pm 1.87$	$0.76\pm 0.65$	$1.49\pm 1.31$	$6.06\pm 5.72$	$90.7$	$95.4$	$37.9$	$84.8$

Table 7. Table 7: Regression results for the SPREAD data trained on the DIR-Lab-4DCT data with elastix using 20, 100, 500 and 2000 iterations. The columns indicate the number of features ( N f subscript 𝑁 𝑓 N_{f} ), the mean absolute error (MAE), the accuracy (Acc) and the F1 score. The sub-indices c, p and w correspond to correct [0,3), poor [3,6) and wrong [6, ∞ \infty ) mm classes, respectively.

	$N_f$	MAE	MAEc	MAEp	MAEw	Acc	F1c	F1p	F1w
Intensity	50	$1.90\pm 3.63$	$1.56\pm 1.49$	$1.26\pm 1.01$	$10.83\pm 14.32$	$71.0$	$82.8$	$21.7$	$48.0$
Registration	108	$1.62\pm 3.59$	$1.23\pm 0.88$	$1.13\pm 0.81$	$11.53\pm 14.60$	$77.1$	$87.4$	$27.7$	$53.9$
Combined	158	$1.73\pm 3.56$	$1.36\pm 0.97$	$1.14\pm 0.83$	$11.30\pm 14.49$	$77.2$	$87.2$	$26.0$	$59.9$

Table 8. Table 8: Regression results for the DIR-Lab-COPDgene data with elastix-advanced and ANTs-BSplineSyN registrations trained on the SPREAD and DIR-Lab-4DCT data. The columns indicate the number of features ( N f subscript 𝑁 𝑓 N_{f} ), the mean absolute error (MAE), the accuracy (Acc) and the F1 score. The sub-indices c, p and w correspond to correct [0,3), poor [3,6) and wrong [6, ∞ \infty ) mm classes, respectively.

	$N_f$	MAE	MAEc	MAEp	MAEw	Acc	F1c	F1p	F1w
elastix-advanced
Intensity	50	$2.17\pm 2.34$	$1.69\pm 1.35$	$2.81\pm 2.66$	$5.15\pm 4.44$	$64.2$	$77.9$	$20.8$	$64.6$
Registration	108	$1.84\pm 2.50$	$1.31\pm 1.66$	$2.22\pm 2.12$	$5.36\pm 4.29$	$76.0$	$87.6$	$29.9$	$57.6$
Combined	158	$1.86\pm 2.05$	$1.50\pm 1.16$	$1.92\pm 1.80$	$4.48\pm 4.21$	$75.3$	$86.9$	$29.5$	$66.1$
ANTs-BSplineSyN
Intensity	50	$2.03\pm 2.01$	$1.80\pm 1.25$	$2.20\pm 2.27$	$5.30\pm 5.38$	$57.3$	$71.6$	$14.2$	$62.7$
Registration	108	$1.71\pm 2.39$	$1.43\pm 1.98$	$2.56\pm 2.01$	$5.06\pm 4.67$	$72.8$	$85.5$	$17.3$	$38.5$
Combined	158	$1.73\pm 1.80$	$1.52\pm 1.23$	$2.22\pm 2.27$	$4.45\pm 4.40$	$76.5$	$87.3$	$20.4$	$59.7$

Equations

$$\displaystyle\widehat{\bm{T}}=\text{arg}\,\underset{\bm{T}}{\text{min}}\,\mathcal{C}\big{(}\bm{T};I_{F},I_{M}\big{)},$$

$$\overline{\mathbf{T}}$$

$$\mathrm{std}\mathbf{T}$$

$$\overline{{\mathbf{T}}^{\mathrm{L}}}$$

$$\mathrm{std}{\mathbf{T}}^{\mathrm{L}}$$

$$\mathcal{E}(\mathbf{T})$$

$$\mathcal{E}({\mathbf{T}}^{\mathrm{L}})$$

$$\mathrm{CV}{\mathrm{H}}^{\mathrm{B}\times \mathrm{B}}$$

$$\mathrm{NMI}$$

$$\mathrm{PMI}$$

$$\displaystyle\mathrm{Imp}(x_{i})=\frac{1}{N_{t}}\sum_{t=1}^{N_{t}}\bigg{(}\underset{j\in\Omega}{\mathrm{MSE}}\Big{(}\hat{y}_{\pi_{i}j},y_{j}\Big{)}-\underset{j\in\Omega}{\mathrm{MSE}}\Big{(}\hat{y_{j}},y_{j}\Big{)}\bigg{)},$$

$$y=\Vert {\mathbf{T}}^{\mathrm{b}}({\mathbf{x}}_F)-{\mathbf{x}}_M{\Vert }_2,$$

$$\mathrm{MAE}$$

Full text

††journal: Medical Image Analysis

1 Introduction

Image registration is the task of finding the optimal spatial transformation between two or more images. In most registration methods, no assessment of the registration quality is provided, and simply the result is returned. Evaluation of the registration is devolved to human experts, which is very time-consuming and prone to inter-observer errors as well as human fatigue [Murphy et al., 2011b]. Automatic quantitative error prediction of registration would decrease quality assessment time and can provide information about the registration uncertainty. Many medical pipelines are based on registered images and it is important to know the uncertainty of registration before continuing to a next phase in order to prevent accumulation of errors. For example, in online adaptive radiotherapy daily contouring of the tumor and organs-at-risk can be performed with the help of image registration [Thörnqvist et al., 2010]. In this task, quality assessment (QA) is mandatory to ensure patient safety. In addition, the accumulation of delivered dose over several treatment fractions is also impacted by the quality of registration [Murphy et al., 2012, Tilly et al., 2013, Veiga et al., 2015]. Registration quality therefore has to be checked before the treatment starts. Visualizing the error of registration can also be directly helpful in medical applications before making a clinical decision. Smit et al. [2017] localized autonomic pelvic nerves by registering a pre-operative MRI scan to an atlas model that includes nerve information. These nerves are not visible in the MRI scans and are prone to be damaged during a surgical procedure. Utilizing registration uncertainty yield better visualization of the autonomic nerves.

Refinement of registration is another important application of automatic error prediction. Muenzing et al. [2014] improved registration by focusing only on regions with high registration error and discarding pixels which are aligned correctly. Registration refinement can also be done with the feedback of human experts by manually adding several corresponding landmarks [Gunay et al., 2017].

Schlachter et al. [2016] did a comprehensive study on visualization of registration quality with the help of three radiation oncologists on the DIR-Lab-COPDgene data, which has a slice thickness of 2.5 mm. The [average, maximum] TRE of the landmarks that were rated to be of acceptable registration quality with the conventional visualization method (checkerboard visualization and color blended) was [2.3, 6.9] mm, while with the best visualization method (histogram intersection) [1.8, 3.3] mm was achieved.

A few methods have been proposed to detect the misalignment of a pair of images with the purpose to refine the registration result. Rohde et al. [2003] proposed to use the gradient of the cost function to detect which region in the image pair is poorly registered and potentially can be improved. Schnabel et al. [2001] suggested to refine the registration by increasing the number of registration parameters in regions with high local entropy, or with high local variation in the intensity or with relatively steep cost function. In another work, analyzing the shape of the cost function around each voxel was used to estimate the confidence of registration [Saygili et al., 2016]. Park et al. [2004] used normalized local mutual information to find poorly aligned regions in order to increase the number of registration parameters. Forsberg et al. [2011] utilized the outer product of the intensity gradient as an uncertainty measure in multi-channel diffeomorphic Demons registration. Although the mentioned metrics can be used to improve the image registration, it has not been shown how these metrics are correlated with the image registration error.

Several methods exploit continuous probabilistic image registration by utilizing Bayesian inference to achieve an intrinsic transformation uncertainty measure [Risholm et al., 2013, Simpson et al., 2015]. However, it has been shown that there is no clear statistical correlation between transformation uncertainty and registration uncertainty [Luo et al., 2017]. The transformation and corresponding label (of a pair of images) are two random variables and it is not possible to quantify the uncertainty of the corresponding label by the summary statistics of the transformation. Another downside of these methods is that they can only be used for the specific paradigm of Bayesian registration.

Some methods are based on the consistency of multiple registrations between a group of images [Datteri and Dawant, 2012, Gass et al., 2015], but these methods cannot be used in pairwise registrations.

In the stochastic approaches, Kybic [2010] suggested to perform multiple registrations with random sampling of pixels with replacement. He found a correlation between the true registration error and the variation of the 2D translational parameters. The method was not extended to 3D and to nonrigid registration. Hub et al. [2009] calculated the local mean square intensity difference multiple times by perturbing the B-spline grid. They showed that the maximum change of the dissimilarity metric in a local region is correlated with the registration error in that region. The drawback of this method is that it is not efficient in homogeneous areas [Hub and Karger, 2013]. In a related work they showed that the variance of the final deformation vector field (DVF) is related to the registration error [Hub and Karger, 2013], using the Demons algorithm. However, to find large misalignment a large search region is needed.

In this paper, we turn our attention to methods capable of learning the registration error allowing to take advantage of multiple features related to registration uncertainty within a single framework. Muenzing et al. [2012] casted the registration assessment task to a classification problem with three categories (wrong, poor and correct registrations). In their method, they mostly utilize intensity-based features, except for the determinant of the Jacobian of the transformation. Although their training samples consist of manually selected landmarks, later they showed that assessing registration in all regions is possible by interpolation [Muenzing et al., 2014].

In our paper, instead of casting the uncertainty estimation task to a classification problem, we formulate it as a regression problem. To the best of our knowledge, in the field of continuous prediction of 3D registration error, Lotfi et al. [2013] only tested their method on artificially deformed images. Recently Eppenhof and Pluim [2017] estimated the registration error by utilizing convolutional neural networks. Only preliminary results were available for synthetic 3D data.

We explore several features related to the uncertainty of the registration transformation as well as related to intensity. All features are calculated in physical units, i.e. mm, which makes the system independent of voxel size. Finally, features are combined by using regression forests. The proposed method is applied and evaluated on chest CT scans. This work is an extension of Sokooti et al. [2016] with updated methodology and substantially extended evaluation.

2 Methods

2.1 System overview

A block diagram of the proposed algorithm is shown in Fig. 1. The system has two inputs: a fixed image $I_{F}$ and a moving image $I_{M}$. Several registration-based and intensity-based features are generated. A regression forests (RF) is then trained from all features to estimate the registration error.

The proposed system is trained to predict residual distances $y$ (registration errors) obtained from a set of semi-automatically established corresponding landmarks. During evaluation, the prediction result $\hat{y}$ is compared with errors obtained from an independent set of ground truth landmarks, using cross-validation. The proposed system therefore estimates registration errors in physical units, i.e. mm. More information about the ground truth is available in Section 3.1. Details of the features are elaborated in Section 2.3.

2.2 Registration

Registration can be formulated as an optimization problem in which the cost function $\mathcal{C}$ is minimized with respect to $\bm{T}$:

[TABLE]

where $\bm{T}$ denotes the transformation. The optimization is usually solved by an iterative method embedded in a multi-resolution setting. A registration can be initialized by an initial transform $\bm{T}^{\mathrm{ini}}$.

2.3 Features and pooling

The features we used in our system, consist of several registration-based as well as intensity-based features. Some features are intrinsically capable to be calculated over differently sized local boxes, for others, a pool of features is created by computing local averages and maxima afterwards. The features used in this paper are listed in Table 1. We propose the following features:

2.3.1 Registration-based features

Variation of deformation vector field ($\operatorname{\mathrm{std}}{\bm{T}}$): The final solution of an iterative optimization problem can be influenced by the initial parameters. If in a region the cost function has multiple local minima or is semi-flat, a slight change in the initial parameters can lead to a different solution. In contrast, in areas where the cost function is well-defined, variations in the initial state are expected to have much less effect on the final solution. A flow chart of the described feature is available in Fig. 2. Given $P$ random initial transformations $\bm{T}^{\mathrm{ini}}_{i}$, $i\in\{1,\ldots,P\},$ that are used as initializations of the registration algorithm from Eq. (1), the variation in the final transformation results $\widehat{\bm{T}}_{i}$ is a surrogate for the precision of the registration. We propose to use the standard deviation $\operatorname{\mathrm{std}}{\bm{T}}$ of those final transformations as a feature:

[TABLE]

In this work, the initial transformations $\bm{T}^{\mathrm{ini}}_{i}$ are created by uniformly distributed offsets in the range $[-2,2]\,$mm to all B-spline coefficients. The offset range is chosen to be relatively small in comparison to the B-spline grid spacing in order to avoid unrealistic deformation. An example of $\operatorname{\mathrm{std}}{\bm{T}}$ in a synthetically deformed image is given in Fig. 3(a).

Instead of perturbing the initial state of the registration, it is also possible to first perform the registration without any manipulated initial state, resulting in a transformation $\bm{T}^{\mathrm{b}}$ [Klein et al., 2009]. Then, random offsets $\bm{T}^{\mathrm{offset}}_{i}$ are added to $\bm{T}^{\mathrm{b}}$ after which another registration is performed, resulting in $\widehat{\bm{T}_{i}^{\mathrm{L}}}$. This is close to the work of Hub and Karger [2013], and approximately measures the concavity of the cost function. The feature $\operatorname{\mathrm{std}}{\bm{T}^{\mathrm{L}}}$ is then derived akin to Eq. (3):

[TABLE]

It is expected that $\operatorname{\mathrm{std}}{\bm{T}^{\mathrm{L}}}$ is small in regions where the cost function is concave, as by adding small offsets $\bm{T}^{\mathrm{offset}}_{i}$ to the parameters, it can still move back to the previous optimal point. A flow chart of $\operatorname{\mathrm{std}}{\bm{T}^{\mathrm{L}}}$ is shown in Fig 2. $\operatorname{\mathrm{std}}{\bm{T}^{\mathrm{L}}}$ is calculated using the same setting as $\operatorname{\mathrm{std}}{\bm{T}}$, except that only one resolution is used.

If the difference between $\overline{\bm{T}}$ and $\bm{T}^{\mathrm{b}}$ is relatively large, regions indicating a small $\operatorname{\mathrm{std}}{\bm{T}}$ are still potentially regions of low registration quality. We then consider the bias $\mathcal{E}(\bm{T})$ and $\mathcal{E}(\bm{T}^{\mathrm{L}})$ as complementary features to $\operatorname{\mathrm{std}}{\bm{T}}$ and $\operatorname{\mathrm{std}}{\bm{T}^{\mathrm{L}}}$ computed by:

[TABLE]

Coefficient of variation of joint histograms (CVH): Multiple registration results can be used to extract additional information from the matched intensity patterns of the images. Given a fixed image $I_{F}$ and a registration sub-result $I_{M}(\bm{T}_{i})$, we calculate their joint histogram $\mathrm{H}_{i},\forall i$. For identical sub-registrations, all resulting joint histograms are equal. Variation in the joint histograms implies registration uncertainty as a surrogate for registration error. The coefficient of variation of the joint histograms is calculated by dividing the standard deviation of all joint histograms over the average, $\overline{\mathrm{H}}$, of them. This normalization is done to compensate for large differences between the elements of $\overline{\mathrm{H}}$. We obtain the CVH in histogram space as follows:

[TABLE]

where B is the number of histogram bins, and $\epsilon$ a constant to avoid division by zero. In the experiments we set $\epsilon$ to 5. The CVHB×B in histogram space is subsequently transferred to the spatial domain, by assigning voxels $x$ with a particular intensity combination \big{(}I_{F}(x),I_{M}(\bm{T}^{\mathrm{b}}(x))\big{)} the corresponding value from CVHB×B, resulting in the final CVH feature with size equal to the fixed image. Note that the CVH can be used in a multi-modality setting, like the previous features. An example of the CVH on a synthetically deformed image is given in Fig. 3(b).

Determinant of the Jacobian (Jac): Jac measures the relative local volume change. This can point to poor registration quality in case of very large ($\mathrm{Jac}\gg 1$) or very small ($\mathrm{Jac}\ll 1$) values, or discontinuous transformations in case of a negative value ($\mathrm{Jac}<0$). In the experiments, the determinant of the Jacobian of $\bm{T}^{\mathrm{b}}$ is used.

2.3.2 Intensity-based features

MIND: The Modality Independent Neighborhood Descriptor (MIND) was introduced by Heinrich et al. [2012] in order to register multi-modal images. In this local self-similarity metric, a patch is considered to compare intensities between fixed and moving images. Finally, the sum of absolute differences between the MIND vector of $I_{F}$ and that of $I_{M}(\bm{T}^{\mathrm{b}})$ is computed. We calculate MIND with a sparse patch including 82 voxels inside a $[7\times 7\times 3]$ box, which is approximately physically isotropic for the data used in the experiments (see Fig. 4).

Local normalized mutual information: Mutual information is used as an entropy-based similarity measure of two images. Similar to Muenzing et al. [2012] we use the following definitions for local normalized mutual information:

[TABLE]

Both metrics are calculated over 8 differently sized boxes: [5, 10, 15, 20, 25, 30, 35, 40] mm. Two strategies for the selection of the number of bins are used, one uses a constant value $\mathrm{B_{C}}$, the other strategy depends on the number of samples $\mathrm{|B|}=log_{2}(n)+1$, in which $n$ is the number of samples in each box. The notations $\mathrm{NMIS}$ and $\mathrm{PMIS}$ indicate mutual information calculated with the latter strategy.

Modality-dependent features: In addition to the modality-independent features from above, we consider the use of several modality-dependent features. In the experiments we assess their contributed value. Similar to Muenzing et al. [2012] the squared intensity difference (SID) and the gradient of intensity difference (GID) are computed using Gaussian (derivative) operators with standard deviations of [0.5, 1, 2, 4, 8, 16] mm. Normalized correlation (NC) is calculated within boxes of size [5, 10, 15, 20, 25, 30, 35, 40] mm akin to Muenzing et al. [2012].

2.3.3 Pooling

In order to reduce discontinuities and improve interaction with other features, the total set of features is increased by generating a pool from those mother features by calculating averages and maxima over them using differently sized boxes. The features MI, SID, GID and NC are inherently computed over differently sized local regions. The features $\mathrm{MIND}$, $\operatorname{\mathrm{std}}\bm{T}$, $\operatorname{\mathrm{std}}\bm{T}^{\mathrm{L}}$, CVH, $\mathcal{E}(\bm{T})$, $\mathcal{E}(\bm{T}^{\mathrm{L}})$ and Jac are calculated in a voxel-based fashion, and then pooled afterwards. Average and maximum pooling is performed with box sizes of $[2,5,10,15,20,25,30,35,40]$ mm. As a result, for each feature we obtain a pool of 18 features: 9 from box averages and 9 from box maxima. The average-pooling is done efficiently by the help of integral images introduced by Viola and Jones [2004]. A list of the proposed mother features together with the number of derived features $N_{f}$ are given in Table 1.

2.4 Regression forests

Random forests were introduced by Breiman [2001] by extending the idea of bagging. The forests consist of several weak learners (trees) which are combined in an efficient fashion. Each tree is started from a node and continues splitting until reaching certain criteria. In contrast to bagging, splitting is performed with a random subset of features which makes the training phase faster and reduces correlation between trees, consequently decreasing the forest error rate. The reason that we chose the random forest is that it can handle data without preprocessing. For instance rescaling of data, outlier removal and selection of features are not necessary in random forests. In addition, random forest are efficient to train and fast at runtime.

Random forests have the capability to calculate the importance of each feature with a little additional computation, which shows the contribution of each feature to the forest. Training of each tree is based on a bootstrap of all samples, and the so-called out-of-bootstrap samples $\Omega$ are used to compute the importance of a feature $x_{i}$. Importance is then defined as the difference between the mean square error (MSE) before and after a permutation of this feature:

[TABLE]

where $y_{j}$ is the real value, $\hat{y_{j}}$ the predicted value from the regression, $\hat{y}_{\pi_{i}j}$ the predicted value when permuting feature $i$, and $N_{t}$ the number of trees.

In this work, random forests are trained with different combinations of the proposed features (see Table 1). The dependent variable $y$ is the registration error in mm, which is described in Section 3.1.

3 Experiments and results

3.1 Materials and ground truth

The SPREAD [Stolk et al., 2007] DIR-Lab-4DCT [Castillo et al., 2009] and DIR-Lab-COPDgene [Castillo et al., 2013] databases have been used in this study. In the SPREAD study, there are 21 pairs of 3D follow-up lung CT images. Each patient in this database has a baseline and a follow-up image (which is taken after 30 months) both in inhale phase. The age of the patients ranges from 49 to 78 years old. The average size of the images is $446\times 315\times 129$ with an average voxel size of $0.78\times 0.78\times 2.50$ mm. In each pair of images, about 100 well-distributed corresponding landmarks were previously selected [Staring et al., 2014] semi-automatically on distinctive locations [Murphy et al., 2011a].

From the DIR-Lab-4DCT data, five cases (4DCT1 to 4DCT5) are selected with each five phases between maximum inhalation and exhalation. The average image size is $256\times 256\times 103$ with an average voxel size of $1.10\times 1.10\times 2.50$ mm. Each scan has 75 corresponding landmarks annotated. Ten cases with severe breathing disorders are available via the DIR-Lab-COPDgene database. The images are taken in inhale and exhale phases. In total, 300 landmarks are annotated. The average image size and the average voxel size are $512\times 512\times 120$ and $0.64\times 0.64\times 2.50$ mm, respectively.

Accuracy of the registration can be defined as the residual Euclidean distance after registration between the corresponding landmarks:

[TABLE]

with $\bm{x}_{F}$ and $\bm{x}_{M}$ the corresponding landmark locations. Based on the idea that the registration error is smooth, we include voxels from a small local neighborhood around the landmarks to increase the total set of available landmarks. In this small neighborhood we assume that the registration error is equal to the error at the center of the neighborhood. This assumption seems reasonable for smooth transformations and within a small region. The neighborhood size is chosen as $10\times 10\times 7.5\;\mathrm{mm}$, which is approximately equivalent to the final grid spacing of the B-spline registration (see Fig. 5).

The core software is written in Python. The feature pooling is performed with a C++ program [Glocker et al., 2014] and the regression forest is calculated with the help of the Scikit-learn package [Pedregosa et al., 2011]. All registrations are performed by elastix [Klein et al., 2010]. Detailed registration setting can be found in the elastix parameter file database (elastix.isi.uu.nl, entry par0049). The code is publicly available via github.com/hsokooti/regun.

3.2 Evaluation measures

In the SPREAD database, we employ 10 cross-validations by randomly splitting the data in 15 image pairs for training and the remaining 6 pairs for testing. To evaluate the regression performance, the mean absolute error (MAE) of the real registration error $y_{i}$ and the estimated one $\hat{y}_{i}$ is calculated over the neighborhood of the landmarks by:

[TABLE]

To further detail the regression performance, the MAE is subdivided into three categories: $\mathrm{MAE_{c}}$, $\mathrm{MAE_{p}}$ and $\mathrm{MAE_{w}}$ with $y$ in $[0,3)$, $[3,6)$ and $[6,\infty)$ mm, corresponding to correct, poor and wrong registration, similar to Muenzing et al. [2012]. We then do the same for $\hat{y}_{i}$, and report the accuracy and F1 score for classifying the registration error in these three categories.

3.3 Parameter selection

The RF is trained using 100 trees with a maximum tree depth of 9, while at least 5 samples remain in the leaf nodes. At each splitting node, $m$ features are randomly selected. We set $m$ to the square root of the total number of features in that experiment, which performed slightly better than $m=\text{(number of features)}/3$ [Liaw et al., 2002]. The total number of registrations $P$ is chosen as 20 to ensure that the estimation of $\operatorname{\mathrm{std}}\bm{T}$ does not change considerably when increasing the number of registrations [Sokooti et al., 2016].

3.4 Reference registration error set

For the SPREAD and the DIR-Lab-4DCT study, registrations are based on free-form deformations by B-splines [Rueckert et al., 1999]. The cost function is mutual information, which is optimized by adaptive stochastic gradient descent. We used three resolutions with a final B-spline grid spacing of $[10,10,10]$ mm. We collect samples by performing four different registrations using 20, 100, 500 and 2000 iterations, respectively. All other registration settings remain the same in these registrations. By varying the number of iterations we increase the variation in the samples, as well as the training size. Table 2 gives the distribution of reference registration errors in each database. As expected, increasing the number of iterations shifts the distribution towards the “correct” registration category. The maximum registration error is 81.8 mm in the SPREAD database, 17.6 mm in the DIR-Lab-4DCT database.

Since the a priori distribution of registration errors is imbalanced, with much more samples in the “correct” category, we perform the following balancing step during training. For landmarks that fall in the category “correct”, we only add samples from a smaller neighborhood of $5\times 5\times 2.5\;\mathrm{mm}$ instead of the $10\times 10\times 7.5\;\mathrm{mm}$ neighborhoods used for landmarks in the categories “poor” and “wrong”. The distribution of reference registration errors of the training samples is shown in Table 3.

For the DIR-Lab-COPDgene study, more advanced settings of the registration are used. In this experiment, samples are taken only on the landmark locations. More details are given in Section 3.5.8. The maximum registration error in this data is 31.5 mm.

3.5 Experiments

3.5.1 Single feature performance in SPREAD

The proposed features are described in Section 2.3 and summarized in Table 1. To investigate the strength of the individual features, we trained the random forest with only a single feature with pooling. By comparing the MAE results in Table 4, it can be seen that MIND, $\operatorname{\mathrm{std}}\bm{T}^{\mathrm{L}}$ and SID$\&$GID are the best single features in the categories Intensity, Registration and Modality-dependent, respectively.

3.5.2 Combined features performance

Instead of using only a single feature, several combinations of features are used to build the RFs:

Intensity: Combination of all modality-independent intensity features: $\mathrm{MIND}$ and MI (50 features).

• 2.

Registration: Combination of all registration features: $\operatorname{\mathrm{std}}\bm{T}$, $\operatorname{\mathrm{std}}\bm{T}^{\mathrm{L}}$, CVH, $\mathcal{E}(\bm{T})$, $\mathcal{E}(\bm{T}^{\mathrm{L}})$ and Jac (108 features).

• 3.

Combined: Combination of both intensity and registration features (158 features).

All results are available in Table 5. By combining features from both the registration and modality-independent intensity category, improvements were obtained in all evaluation measures.

The result of the regression with combined features is detailed in Fig. 6(a), which shows the real error (solid blue line) against the predicted error, sorted from small to large. In Fig. 6(b) we grouped the real errors in the three categories, each category showing a box-plot of the predicted errors. Intuitively, a smaller overlap between the boxes represents a better regression.

3.5.3 Including modality-dependent features

We consider adding the combination of three modality-dependent features to the combined feature set (Combined+MD): NC, SID and GID. In both databases, if we add the modality-dependent features (see Table 5), negligible differences are observed. Therefore, to keep the feature set small and modality-independent, we select the “combined features” class without the modality-dependent features as the final system in the remainder of this paper.

3.5.4 The effect of pooling

To examine the effect of pooling, we perform an experiment without pooling on the combined feature set. We only calculate PMIS within a box size of 15 mm in this experiment. From Table 5 the benefit of pooling can be observed.

3.5.5 Alternative regression methods

In this section, we compare RF regression with linear regression (LR) and neural networks (NN). Feature normalization is done for both regressors. We utilized neural networks with three hidden layers of 1024, 512 and 256 units each. ReLU is used as an activation function and Huber is utilized as a loss function. Table 5 gives the results of these experiments. The performance of neural networks is on par with random forests. However, the results of linear regression are not comparable to that of random forests, both in MAE and accuracy.

3.5.6 Feature importance

The feature importance, see Eq. (9), is displayed in Fig. 7. It shows that MIND and MI are the features contributing most to the RF performance, followed by $\operatorname{\mathrm{std}}{\bm{T}}$, $\operatorname{\mathrm{std}}\bm{T}^{\mathrm{L}}$ and CVH.

The feature importance using a different number of iterations is shown in Fig. 8. The contribution of all intensity features stay the same in all experiments, while some of the registration features contribute differently with respect to the number of iterations. For instance, the importance of $\operatorname{\mathrm{std}}\bm{T}$ and CVH increase with increasing the number of iterations. The features $\operatorname{\mathrm{std}}\bm{T}^{\mathrm{L}}$ and $\mathcal{E}(\bm{T}^{\mathrm{L}})$ play important roles when the number of iterations is not enough for registration convergence.

3.5.7 Excluding a single feature

To further investigate the importance of the several features, we additionally perform an experiment where we leave one feature out of the combined feature set. The results are reported in Table 6. In these experiments, feature redundancy can be found. For instance, MI has a large importance values in random forests, but if we leave that feature out, other features can compensate for that.

3.5.8 Inter-database validation

To study the generalizability of the proposed system, instead of cross-validation on a single database, we perform training on the DIR-Lab-4DCT database and test it on the SPREAD database. As mentioned before, the SPREAD database consists of only inhale images but the DIR-Lab-4DCT database has images from inhale to exhale phases. Therefore, this makes the DIR-Lab-4DCT more suitable for training. The result of this experiment is available in Table 7. Once more, we can draw the conclusion that by combining both intensity and registration-based features, the regression performance can be improved. In contrast to the SPREAD experiment, this time it is observed that the registration features perform better than the intensity features.

To further evaluate the generalizability of the proposed method, we test it for different registration methods on a third independent test set, the DIR-Lab-COPDgene dataset. The regression forest is trained on a combination of the SPREAD and DIR-Lab-4DCT data. We evaluate two registration algorithms that achieved excellent performance in the EMPIRE10 challenge [Murphy et al., 2011b], i.e. the ANTs registration package [Avants et al., 2009, Tustison et al., 2013] and elastix with advanced settings [Staring et al., 2010].

Prior to deformable registration we perform an affine registration using 5 resolutions and utilizing torso masks. For the deformable registration we use settings similar to the ones used in the EMPIRE10 challenge, specifically:

ANTs-BSplineSyN: With respect to the EMPIRE10 challenge we increased the number of iterations to 1000 for each of the 4 resolutions, using a 10% sampling rate. This improved the performance on our data and considerably reduced the calculation time. As suggested in Tustison et al. [2013], several preprocessing steps are used, including masking out the lungs, and inverting the image intensities and rescaling them between 0 and 1. Further settings include: registration model: symmetric diffeomorphic; dissimilarity metric: local cross correlation; number of resolutions: 4; maximum number of iterations: 1000; sampling: 10% random samples; convergence threshold: 1e-6. The average TRE on DIR-Lab-COPDgene is $1.90\pm 2.86$ mm.

elastix-advanced: Settings are adopted from Staring et al. [2010]. The most important ones are: registration model: B-spline; dissimilarity metric: normalized correlation; number of resolutions: 6; number of iterations: 1000; sampling: 2000 random samples; B-spline grid spacing: [5, 5, 5] mm. The average TRE with this setting is $3.39\pm 4.30$ mm on the DIR-Lab-COPDgene dataset.

Detailed parameter files for both registration methods are available via elastix.isi.uu.nl (entry par0049) and github.com/hsokooti/regun. The calculation time of ANTs was about 60 hours per registration, comparing to 12 minutes for elastix.

In this experiment, the evaluation is performed only on the landmarks locations, where Table 2 displays the distribution of reference registration errors during testing. The results of the experiments are given in Table 8. A scatter plot is also depicted in Fig. 9. Similar to the previous inter-database experiment (Table 7), the MAE and accuracy of the registration features are slightly better than the MAE and accuracy of the intensity-based features. However, intensity features obtained better classification score in the wrong category. We conclude that the proposed method indeed generalizes to different settings of the same method (elastix-advanced), as well as registration methods with quite a different underlying transformation model (ANTs-BSplineSyN, which uses a symmetric diffeomorphic model).

4 Discussion

A system for quantitative error prediction of medical image registration is proposed and it is quantitatively evaluated on multiple chest CT datasets.

4.1 Features

In the intra-database (SPREAD) validation, it is observed that the single MIND feature can perform almost as good as the overall combined system. By adding MI and registration features, the results slightly improved. Muenzing et al. [2012] did not consider MIND in their feature set and found that the most important single features in their classification experiments are mutual information and Gaussian intensity, whereas, based on Table 4 these features are less important than MIND in our experiments. Furthermore, the calculation time of MI for the whole image is about 3 h, as opposed to the calculation time of MIND, which is about 8 min ($\sim$3 min MIND + $\sim$5 min pooling). Although less accurate, it is possible to reduce the calculation time of the MI feature by calculating it over a single window and then aggregate by pooling.

The modality-dependent intensity features do not increase regression accuracy on the data used in our paper. Consequently more generally applicable modality-independent features can be used, even for mono-modal problems.

Table 5, 7 and 8 together suggests that features in the intensity and registration categories provide complementary information, and that a better system can be obtained in terms of MAE and accuracy by considering both intensity and registration-derived features.

The intensity features were better predictors than the registration features in the intra-database experiment. However, in the inter-database experiment, the registration features outperform intensity features in terms of total accuracy and MAE. The same observation can be made for the average F1 score in the inter-database experiments using elastix (See Table 7, 8). For ANTs (Table 8), the average F1 score of the intensity-based features was slightly higher than that of the registration-based features.

The registration features contribute differently with respect to the number of iterations (See Fig. 8). The features $\operatorname{\mathrm{std}}\bm{T}^{\mathrm{L}}$ and $\mathcal{E}(\bm{T}^{\mathrm{L}})$ play important roles when the number of iterations is not enough for convergence. When the number of iterations increases, the contribution of $\operatorname{\mathrm{std}}\bm{T}$ and CVH go up. In the work of Muenzing et al. [2012], only one registration feature, Jac, was used and they reported that the impact of this feature is relatively low in comparison with intensity-based features. We observed the same result for Jac, but it should be pointed out that the range of Jac in our database was [0.3, 3.9] so voxels with negative or very large values were not encountered.

Feature pooling improves the regression results in all evaluation measures, due to the addition of contextual information. In some features like $\operatorname{\mathrm{std}}\bm{T}$, average pooling contributed more to the regression performance, while in features like CVH, maximum pooling had a higher importance value (See Fig. 8(d)).

As can be seen in Table 6, the proposed combined system has redundant features. Hence, by removing a single feature, the system is still able to predict the registration error with almost equal MAE as the total system. However, removing these features may decrease the generalizability of the system. For example, looking at the feature $\mathcal{E}(\bm{T}^{\mathrm{L}})$ in Fig. 7 we see that its contribution is relatively small overall. However, in Fig. 8 it can be seen that while it is less important for better registration results (100, 500 and 2000 iterations), it is still important for poor registration results (20 iterations).

Considering the results in both intra and inter-database experiments (Table 5, 7 and 8), the conclusion to be drawn is that the proposed combined feature sets is general and robust.

4.2 Quantitative validation

Commonly, in image registration tasks, the distribution of registration errors is not balanced as can be seen in Table 3.

In the SPREAD experiment, Table 5 reports that in the combined experiment, the MAE of the correct and poor classes are $0.76\pm 0.65$ and $1.47\pm 1.22$, respectively. On the contrary, the MAE of the wrong class is $6.12\pm 5.64$. It is expected that the regression error of values of the wrong class is relatively larger than that of the other classes. However, it should be emphasized that only $3.9\%$ of samples are available to make a regression model between 6 and 81.8 mm. We tried to add more samples and make the distribution more balanced by performing registrations with different number of iterations, but there is still room for improvement for the wrong class by adding more samples and data.

In terms of classification, we obtained F1 scores of $95.4\%$, $38.1\%$ and $84.4\%$ in the classes correct, poor and wrong, respectively (Table 5). For the wrong class, which is arguably most important for clinical application, the precision and recall are $84.6\%$ and $84.3\%$, respectively. This means that $84.6\%$ of all samples predicted to be over 6 mm are correct and the proposed method caught $84.3\%$ of larger registration errors. Muenzing et al. [2012] obtained F1 scores of $95.3\%$, $73.8\%$ and $86.6\%$ in the classes [0,2), [2,5) and [5, $\infty$) mm. They achieved a better F1 score in the poor class and they also reported zero overlap between the correct and wrong classes. However, the comparison between the two methods is not easy because of the differences in the data. For example, the slice thickness in SPREAD is 2.5 mm, while it is 1 mm for Muenzing’s data, which may affect the performance especially in the poor class. Moreover, we generated the classes by thresholding the regression values. Thus, the forests are optimized for regression not for classification.

4.3 Qualitative validation

Muenzing et al. [2014] generated an uncertainty map by spatial interpolation of landmark-based quality estimates. On the contrary, our proposed system, which is trained on landmark locations, can be applied in all regions of the image. We showed this for two example images, see Fig. 5. It can be easily visualized that in the blue region, images are matched correctly. On the other hand, by tracking the vessels in the red region misalignment can be seen. Another note about the prediction is that there are no abrupt changes, and error varies smoothly from blue to yellow and then red, even though the error is predicted for each voxel independently.

Another example is given in Fig. 10(a-d). Although all landmarks indicate that the registration error is small in this slice, the quantitative results found several misregistered regions, which implies that few landmarks may not be sufficient to assess the registration quality of the whole image. In Fig. 10(e, f), it can be observed that the performance in the homogeneous area (left side of the images) is as good as the performance in the area with structure. The main reason of acceptable performance in the homogeneous area is that the training samples consist of landmarks as well as their neighborhood region, which can be homogeneous. Thus, the system is trained both for homogeneous regions and regions with structure.

Another example is given in Fig 10(g, h), where the proposed system is not able to predict the registration error correctly because of a shift in the slice direction.

4.4 Limitations

Discrete optimization: If the optimization method is less or not dependent on the initial state, for instance for discrete optimization methods [Glocker et al., 2008, Heinrich et al., 2016], many of the proposed registration features, which are generated by varying the initial transformation of the registration, are not informative anymore. In such cases, instead of $\operatorname{\mathrm{std}}{\bm{T}}$ or $\operatorname{\mathrm{std}}{\bm{T}^{\mathrm{L}}}$, other measures can be used. For example, by utilizing the adaptive mean-shift algorithm, the local standard deviation of the displacement distribution can be calculated [Heinrich et al., 2016].

Anatomical changes: The proposed method is trained in such a way that any dissimilarity between the fixed and moving images is counted as misalignment in registration. In case of anatomical changes this assumption may be invalid, but typically prior knowledge of the underlying anatomy is required to determine which regions are allowed to be ”misaligned” because of anatomical changes and which are not [Muenzing et al., 2018]. The proposed method highlights all changes, coming from misalignment or from anatomical change.

4.5 Future work

In the proposed method we predict the misalignment as an Euclidean distance in millimeters, rather than a 3D vector representing residual displacement. This is mostly because the features used in the system are not direction-wise, especially the local intensity features. The use of features that include directional information may help the system to be used in predicting the registration error in each direction, which is then effectively a new registration method.

The proposed method was tested on chest CT scans. Since the proposed features are generic and modality-independent, the overall method can in principle be applied to other modality data from other anatomical regions. The performance in such cases however remains to be investigated.

The uncertainty of affine registration is not measured in this work. Defining a gold standard for this mid-phase result is a complex task. However, extending the experiments to other databases where only affine transformations are applicable can be done in the future.

Instead of manually defined features, it is possible to use convolutional neural networks, which can learn features automatically. Eppenhof and Pluim [2017] predicted the Euclidean distance of registration error. Our own work on CNNs for registration [Sokooti et al., 2017] can also be modified to predict registration uncertainty in a direction-wise manner. Both methods are trained only based on intensity, where the current paper shows that registration-derived information still contributes to a better regression. Thus, adding registration information to the neural networks should probably be considered as well.

A larger set of corresponding points annotated more densely throughout the scan could potentially also benefit training of the regression forest. In addition, experimenting on multi-modality data and investigating the contribution of all introduced features on them are future plans of this work.

Finally, the uncertainty map produced by the proposed method may be exploited to improve local registration results.

5 Conclusion

In this paper we proposed a method based on random regression forests to predict registration accuracy on chest CT scans from registration-based as well as intensity-based features. We introduced the variation in registration result from differences in initialization ($\operatorname{\mathrm{std}}{\bm{T}}$) and CVH, which showed high feature importance in several experiments. Registration-based features provided additional information on registration error with respect to intensity-based features.

The regression method was evaluated on data from the SPREAD study and predicted the registration error with a mean absolute error of 1.07 $\pm$ 1.86 mm. The proposed method gained promising results on inter-database validation with a regression error of 1.76 $\pm$ 2.59 mm.

Acknowledgments

This work is financed by the Netherlands Organization for Scientific Research (NWO), project 13351. Ben Glocker received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 757173, project MIRA, ERC-2017-STG). Dr. M.E. Bakker and J. Stolk are acknowledged for providing a ground truth for the SPREAD study data used in this paper. We would like to thank Dr. R. Castillo and T. Guerrero for providing the DIR-Lab database.

6 References