Keywords

1 Introduction

Medical image registration is a vital tool in defining spatial relations between multiple medical images. In medical image registration, performance and outcome of a registration are heavily determined by several user-defined parameters [9] such as penalty weights, number of iterations, smoothing parameters, etc. Optimal values of these parameters are heavily affected by the anatomical region of interest, image acquisition settings, the expected severity of deformations, and the clinical requirements. Therefore, in order to get optimal registration performance for a given image pair, methods providing insight into the effect of these parameters are needed. An exhaustive search over the parameter space could be a solution for this problem when computation time is not an issue, or when only a single parameter needs to be tuned. However, it is not suitable for some clinical applications where the time dedicated to the registration is limited such as image-guided interventions. Moreover, when multiple parameters need to be tuned, an exhaustive search over the parameter space becomes infeasible due to exponentially increasing computation time (curse of dimensionality).

In this study, we introduce a method for efficient parameter tuning based on polynomial chaos expansion (PCE). PCE is used for uncertainty quantification in a variety of fields such as solid mechanics, fluid mechanics, and econometrics. PCE is a stochastic approach to relate a model’s output variability (in this study, registration output) with the variability of a random input (registration tuning parameters) using orthogonal polynomials in a compact representation [2]. It requires sampling of the model output for a modest number of input values, derived from the distributions of the input variables. Using the resulting compact polynomial representation, a user is able to acquire outputs analytically for corresponding user-defined inputs without running any further registration. The sampling of model outputs is done with a few registration runs whose number is considerably lower than those of exhaustive search. The number of sampling points in PCE is much lower than the number of registrations required by an exhaustive search procedure, especially when exploiting sparse grid sampling methods for PCE [5]. Even when multiple inputs are used, the construction of the PCE only requires a limited number of model evaluations [1, 4, 5], making it applicable for time-critical applications. After establishing the PCE model, registration results can be acquired for given input combinations and the user can simulate registration results in real-time. The method is tested on 8 datasets and results indicate that the method is promising for the mentioned use.

2 Method

2.1 Theory

PCE was first introduced by Wiener [10] who showed that series of polynomials of centered normalized Gaussian variables can approximate any second order random variable:

$$\begin{aligned} \mathbf {Y} _{PCE} = \sum _{k=1}^{K} c_k\mathbf {\psi }_k(\mathbf X ), \end{aligned}$$
(1)

where \(\mathbf X \) is a Gaussian random variable vector, \(\mathbf {\psi }_k\) are multivariate Hermite polynomials, \(c_k\) are their coefficients and K defines polynomial order. The \(\mathbf {\psi }_k\)’s with different indexes are orthogonal to each other with respect to the Gaussian measure. Generalization of the PCE to other, non-Gaussian, random variable distributions was presented by Xiu et al. [11]. Each distribution type is associated with a specific polynomial type, as shown for some well-known examples in Table 1. The main task for model construction is to find coefficients \(c_k, k=1\ldots K\) such that \(\mathbf Y _{PCE}\) approximates the true output of the model. In order to construct the PCE model, three steps are followed. First, the distributions of the input variables are determined and corresponding polynomial types are selected. Next, the polynomial order K and the so-called “grid level” (to control the sparseness of the sampling grid [5]) are determined. Subsequently, the sampling of the real function is triggered and the PCE model is constructed by calculating the coefficients \(c_k\). Once the PCE model is constructed, the output can be evaluated for any desired input parameter setting [6].

Table 1. Some random variable distributions and their corresponding polynomials.
Full size table

2.2 Distribution Assumptions for User-Defined Input Parameters

In image registration, penalty terms play an important role in improving registration accuracy. Basically, they are weighted against a (dis)similarity function via their weighting coefficients. As an example, we consider a parametric intensity-based registration method involving two penalty terms, mathematically represented by the following optimization problem:

$$\begin{aligned} \mathbf {\hat{{\mu }}} = \underset{{\mathbf {\mu }}}{\arg \min }{\ }C({\mathbf {\mu }};F,M)+{\omega }_1P_1({\mathbf {\mu }})+{\omega }_2P_2({\mathbf {\mu }}), \end{aligned}$$
(2)

where \(\mathbf {\mu }\) is a vector of transformation parameters, C is an intensity-based dissimilarity measure, F is the fixed image and M is the moving image. \({\omega }_1\) and \({\omega }_2\) are weighting terms for the penalty terms \(P_1\) and \(P_2\), respectively. In the experiments, we chose transform bending energy [8] and point-to-surface [3] penalty terms for \(P_1\) and \(P_2\), respectively. To enable efficient interactive tuning of the parameters \(\omega _1\) and \(\omega _2\), we propose to construct a PCE model that predicts \(\mu \) for every combination of input parameters \(\omega _1\) and \(\omega _2\). Since penalty weighting parameters are usually considered on a logarithmic scale, we introduce the following reparameterization:

$$\begin{aligned} {\omega }_1(x_1) = 2^{6+2x_1},{\omega }_2(x_2)=2^{-9+x_2}, \end{aligned}$$
(3)

In this reparameterization, we were motivated by our prior knowledge on the typical range of these tuning parameters and we model \(x_1\) and \(x_2\) as random variables with a zero mean Gaussian distribution with \(\sigma =2\). Finally, our PCE model can be written as:

$$\begin{aligned} \mathbf {{\mu }}_{PCE} = \sum _{k=1}^{K} c_k\mathbf {\psi }_k(x_1,x_2), \end{aligned}$$
(4)

where \(\mu _{PCE}\) are transform parameters predicted by this PCE model.

2.3 PCE Model Construction

Since PCE is modeled using orthogonal polynomial bases, the PCE coefficients \(c_k\) in (1) can be determined by spectral projections as follows:

$$\begin{aligned} c_k = \dfrac{\langle \mathbf Y (\mathbf X ), \mathbf {\psi }_k(\mathbf X )\rangle }{\langle \mathbf {\psi }_k(\mathbf X ),\mathbf {\psi }_k(\mathbf X ) \rangle }=\dfrac{\int ...\int \mathbf Y (\mathbf X )\mathbf {\psi }_k(\mathbf X )p_\mathbf{x }(\mathbf X )dx_1...dx_N}{\langle \mathbf {\psi }_k(\mathbf X ),\mathbf {\psi }_k(\mathbf X ) \rangle } \end{aligned}$$
(5)

In PCE, the multidimensional integration for spectral projection in the equation (5) is computationally demanding and poses a problem for the feasibility of PCE in multiple input-output problems. However, the severity of the problem is considerably alleviated by the introduction of a sparse grid based numerical integration technique [5]. The accuracy of PCE is heavily dependent on polynomial order and sparse grid level selection. It increases with the increasing order and level. However, higher grid level also mean more sampling and this causes a trade-off. After selection of feasible polynomial order and grid level, the sampling of model outputs is done at model inputs determined by the input distribution model, polynomial order, and sparse grid level.

In parametric image registration, results are represented by transformation parameters which define how the deformation of images are evolving. In rigid registration, there are 6 transformation parameters; 3 for translation and 3 for rotation. However, for nonrigid registration, the number of transformation parameters representing the deformation field may be in the order of millions. In our experiments, we use a B-spline transformation model [7]. Each transformation parameter is treated as an independent model output and a PCE model specific to that transformation parameter is constructed. After constructing all PCE models, they can be used to instantly obtain a deformation field for any given value of the input tuning parameters (\({\omega }_1\) and \({\omega }_2\) in Eq. 2).

3 Experiments and Results

3.1 Dataset

In the experiments, we used 8 pairs of abdominal CT images taken from 8 patients for ablation interventions. All images were anonymized prior to processing and acquired at Erasmus MC in 2014. Image resolutions vary from 0.71 mm x 0.71 mm to 0.84 mm x 0.84 mm (501 \(\times \) 492 to 512 \(\times \) 512 in-plane resolution), 3–5 mm slice spacing and 1–2 mm slice thickness (16 to 70 slices). The point-to-surface penalty term requires a set of manual landmarks in the fixed image and a segmentation of the liver in the moving image, as explained in [3], which were available to us.

3.2 Evaluation

For the PCE model, we chose a polynomial order of 4 and a grid sparsity level of 3. With these settings, 21 sampling points (or registrations) at different values of \(\omega _1\) and \(\omega _2\) were required by the PCE. We examined the method for the values of the surrogate variables \({x}_1\) and \({x}_2\) in the range of [-3, 3] with 1 increase at each step. This resulted in 7\(\,\times \,\)7 grid. Therefore, covered ranges of \({\omega }_1\) and \({\omega }_2\) are \(2^0...2^{12}\) and \(2^{-12}...2^{-6}\), respectively. PCE results were examined both qualitatively and quantitatively for the given ranges of \({\omega }_1\) and \({\omega }_2\).

Fig. 1.
figure 1

An axial slice of a moving image after being deformed with PCE-generated (A and C) and real (B and D) transformation parameters for (\(\omega _1 = 2^{8}, \omega _2=2^{-12}\)) (A and B) and (\(\omega _1 = 2^{10}, \omega _2 = 2^{-11}\)) (C and D), respectively.

Full size image
Fig. 2.
figure 2

Transform parameter differences between two consecutive values of \(\omega _1\) (Eq. 7) (A), two consecutive values of \(\omega _2\) (Eq. 8) (B) and approximation error of PCE and real registration results (Eq. 6) (C).

Full size image

3.3 Results

In Fig. 1, we see four panels with PCE-generated (A and C) and original transformation parameters (B and D). Setting 1 (\(\omega _1 = 2^{8}, \omega _2 = 2^{-12}\)) is used to obtain panels A and B, and setting 2 (\(\omega _1 = 2^{10}, \omega _2 = 2^{-11}\)) is used to obtain C and D. Panels A and B show that the transformation parameters predicted by PCE are similar to those obtained by a real registration for setting 1. Panels C and D show the same for setting 2. Comparing panels A/B with C/D, it is clear that settings 1 and 2 resulted in a very different registration outcome. The PCE model successfully imitated the real registration’s response to changes of the tuning parameters \(\omega _1\) and \(\omega _2\).

To quantitatively evaluate the accuracy of the PCE, we computed the error between the PCE-generated (\(\mathbf {\mu }_{pce}(\omega _1,\omega _2)\)) and real (\(\mathbf {\mu }_{re}(\omega _1,\omega _2)\)) transformation parameters:

$$\begin{aligned} {\epsilon }_\mathbf {\mu }(\omega _1,\omega _2) = \dfrac{1}{N}\sum _{j = 1}^{N} |\mu _{pce}^j(\omega _1,\omega _2)-\mu _{re}^j(\omega _1,\omega _2)|, \end{aligned}$$
(6)

where N is number of transformation parameters. In order to quantify the sensitivity of the real registration’s outcome to small changes in \(\omega _1\) and \(\omega _2\), we examined differences between transformation parameters for two consecutive values of \(\omega _1\) and \(\omega _2\), respectively, using the equations below:

$$\begin{aligned} {d}_{\mathbf {\mu }_1}(\omega _1,\omega _2) = \dfrac{1}{N}\sum _{j = 1}^{N} |\mu _{re}^j(\omega _1(x_1+1),\omega _2(x_2))-\mu _{re}^j(\omega _1(x_1),\omega _2(x_2))|, \end{aligned}$$
(7)
$$\begin{aligned} {d}_{\mathbf {\mu }_2}(\omega _1,\omega _2) = \dfrac{1}{N}\sum _{j = 1}^{N} |\mu _{re}^j(\omega _1(x_1),\omega _2(x_2+1))-\mu _{re}^j(\omega _1(x_1),\omega _2(x_2))|. \end{aligned}$$
(8)

Figure 2, presents the error between the real and PCE-generated transformation parameters (C) and discrepancies between values of transform parameters when weighting coefficients of transform bending energy (panel A (see Eq. 7)) and point-to-surface (panel B (see Eq. 8)) penalty terms change. On the panels A and B, high intensity pixels exhibit input values where the real registration is sensitive to changes. Looking at the panel C, it is clear that the PCE method approximates the real registration well even for the input values where the registration is sensitive. Therefore, it can be stated that the PCE model generates quite accurate results for broad ranges of \(\omega _1\) and \(\omega _2\) values with the maximum error (on panel C) of 0.63 mm.

Here, it should be noted that the results presented are averaged over 8 image pairs from different patients.

4 Discussion

In medical image registration, there is a multitude of user-defined parameters which need to be tuned for optimal results. Basically, tuning of those parameters could be done by an exhaustive search of the parameter space, but this would be computationally very demanding, and infeasible in time restricted applications. The problems becomes even worse when tuning multiple parameters, due to the curse of dimensionality. To address this problem, we proposed a PCE based method that constructs an approximation of the registration algorithm from a moderate number of registration executions for different values of input parameters. Based on these registrations, a PCE model is constructed and saved in computer memory to instantly obtain any desired transformation parameters for specific values of input parameters.

The PCE model then allows users to instantly see results for selected input parameter values. A PCE-based visualization tool could be developed allowing users to immediately inspect the deformation results for different parameter values, e.g. by simply dragging a slider, thus enabling efficient and interactive parameter tuning. Moreover, if a quantitative measure of registration accuracy can be computed (e.g., Dice overlap between segmented tissues), the parameters could be optimized automatically using the PCE model. A prerequisite for both approaches is that the PCE model provides a good approximation of the true registration result. Verifying this requirement was exactly the purpose of the experiment presented in Sect. 3

The PCE method may have methodological similarities with usual interpolation methods. However, there are two major differences which may affect accuracy and computational demand drastically. First, in interpolation, sampling is done on a uniform base and without harnessing the distribution of input parameters. Thus, the sampling is done sub-optimally. Second, for multiple inputs, sampling is done on a uniform grid and therefore computational demands may explode due the curse of dimensionality. In contrast, the PCE method determines sampling points on a sparse grid by exploiting prior knowledge on the distribution of inputs. In this way, the computational expenses are reduced considerably.

In the experiments, we evaluated the method on 8 datasets by examining differences between the PCE-generated and real transformation parameters for the same weighting coefficients of two penalty terms in order to show how accurate the method generates transformation parameters. It was shown that the method generates rather accurate transformation parameters compared with the real ones for a broad range of values of input parameters. The maximum discrepancy was 0.63 mm. Further, we also presented how changes in those weighting coefficients affect transformation parameters of real registration so as to find where the registration is sensitive to the changes and whether the PCE method is capable of generating accurate results around sensitive regions. The errors caused by the method were low and this showed that the method can successfully address changes in values of input parameters.

5 Conclusion

In this study, we proposed a polynomial chaos expansion (PCE) based method which can be used to examine and tune user-defined input parameters of a registration in a fast manner. The method generates accurate results and is promising in terms of both timing and performance.