A generalized framework unifying image registration and respiratory motion models and incorporating image reconstruction, for partial image data or full images

To build a respiratory motion model image data are simultaneously acquired with one or more surrogate signals (figure 1). The image data consist of N_i dynamic images, ${{\mathbf{I}}_{1}}\cdots {{\mathbf{I}}_{{{N}_{i}}}}$, each acquired at a different time-point, t, and representing a different respiratory state. These should cover at least one respiratory cycle, although sometimes several cycles will be imaged so that inter-cycle variations can be sampled. A reference-state image, I₀, is also required. This could be one of the dynamic images or it could be another image, e.g. a high quality breath-hold image.

Following the typical approach to building motion models shown in figure 1, the respiratory motion is first determined form the images, and then a correspondence model it fitted relating the motion and the surrogate signals. To determine the motion I₀ is registered to each of the dynamic images, I_t. Each registration is usually performed independently, and results in a spatial transformation parameterized by the motion parameters, $~{{\mathbf{M}}_{t}}=\left[{{m}_{t,1}},\cdots,{{m}_{t,{{N}_{\text{m}}}}}\right]$, where N_m is the number of motion parameters. e.g. M_t could be the parameters of an affine transform, a B-spline control point grid, or a voxel-wise deformation field.

The correspondence model relates M_t to the surrogate signals, ${{\mathbf{S}}_{t}}=\left[{{s}_{t,1}},\cdots,{{s}_{t,{{N}_{\text{s}}}}}\right]$, where N_s is the number of surrogate signals measured at each time-point, t. The model is parameterized by the model parameters, $\mathbf{R}\,={{\mathbf{R}}_{1}}\cdots {{\mathbf{R}}_{{{N}_{r}}}}$, where N_r is the number of model parameters per motion parameter, and ${{\mathbf{R}}_{n}}=\left[{{r}_{1,n}},\cdots,{{r}_{{{N}_{\text{m}}},n}}\right]$, i.e. the total number of model parameters is N_r  ×  N_m. The correspondence model can be represented as a function, F, of S_t and R that results in an estimate of the motion, ${{\mathbf{M}}_{{{F}_{t}}}}$:

$$\begin{eqnarray}&&{{\mathbf{M}}_{{{F}_{t}}}}=F\left({{\mathbf{S}}_{t}},\boldsymbol{}~\mathbf{R}\right)\end{eqnarray}$$

The correspondence model is fitted using a standard technique such as ordinary least squares, which minimizes the squared difference between the motion estimated by the model and the registration results over all time-points:

$$\begin{eqnarray}&&\underset{\mathbf{R}}{{\min}}\,\underset{t}{\sum}\,\|{{\mathbf{M}}_{t}}-{{\mathbf{M}}_{{{F}_{t}}}}\|_{2}^{2}\end{eqnarray}$$

The framework presented in this paper does not follow the typical approach described above. Instead, it directly optimizes the correspondence model parameters, R, on all the dynamic images simultaneously, such that the motion calculated using the correspondence model transforms I₀ to best match the dynamic image data, i.e. ${{\mathbf{M}}_{t}}={{\mathbf{M}}_{{{F}_{t}}}}$.

Most image registration algorithms try to find the optimal values of M_t by minimising a cost function, C_t, which consists of a (dis)similarity term, Sim, and an (optional) constraint term, Con.

$$\begin{eqnarray}&&{{C}_{t}}\left({{\mathbf{I}}_{t}},{{\mathbf{I}}_{0}},{{\mathbf{M}}_{t}}\right)=Sim\left({{\mathbf{I}}_{t}},T\left({{\mathbf{I}}_{0}},{{\mathbf{M}}_{t}}\right)\right)+Con\left({{\mathbf{M}}_{t}}\right)\end{eqnarray}$$

where I_t is the fixed (target) image, I₀ the moving (source) image, and T the function that applies the transformation parameterized by M_t to image I₀. Note, I₀ is used as the moving image here as this is more convenient when using partial image data (see section 2.3). A common approach for optimizing the values of M_t is to estimate the gradient of C_t with respect to M_t, i.e. $\frac{\partial {{C}_{t}}}{\partial {{\mathbf{M}}_{t}}}$, and then use a gradient based optimization such as gradient descent or conjugate gradient.

To directly optimize the correspondence model parameters on the image data the gradient of the cost function with respect to the model parameters, $\frac{\partial {{C}_{t}}}{\partial \mathbf{R}}$, is required. If the correspondence model can be differentiated to give the gradient of the motion parameters with respect to the model parameters, $\frac{\partial {{\mathbf{M}}_{t}}}{\partial \mathbf{R}}$, the chain rule can be applied, i.e.:

$$\begin{eqnarray}&&\frac{\partial {{C}_{t}}}{\partial \mathbf{R}}=\frac{\partial {{\mathbf{M}}_{t}}}{\partial \mathbf{R}}\frac{\partial {{C}_{t}}}{\partial {{\mathbf{M}}_{t}}}\end{eqnarray}$$

Therefore, any correspondence model can be used in this framework as long as it is possible to calculate $\frac{\partial {{\mathbf{M}}_{t}}}{\partial \mathbf{R}}$. This includes all of the popular correspondence models that have been used in the literature, such as linear, polynomial, and B-spline models (McClelland et al 2013). e.g. for a 2nd order polynomial model with one surrogate signal:

$$\begin{eqnarray}&&{{\mathbf{M}}_{t}}={{\mathbf{R}}_{1}}s_{t}^{2}+{{\mathbf{R}}_{2}}{{s}_{t}}+{{\mathbf{R}}_{3}}\end{eqnarray}$$

$$\begin{eqnarray}&&\frac{\partial {{\mathbf{M}}_{t}}}{\partial {{\mathbf{R}}_{1}}}=s_{t}^{2},\boldsymbol{~}\frac{\partial {{\mathbf{M}}_{t}}}{\partial {{\mathbf{R}}_{2}}}={{s}_{t}},\boldsymbol{~}\frac{\partial {{\mathbf{M}}_{t}}}{\partial {{\mathbf{R}}_{3}}}=1\end{eqnarray}$$

The correspondence model parameters are optimized on all of the image data by summing the cost function and gradient from each of the individual images:

$$\begin{eqnarray}&&{{C}_{\text{total}}}=\underset{t=1}{\overset{{{N}_{i}}}{\sum}}\,{{C}_{t}},\quad \frac{\partial {{C}_{\text{total}}}}{\partial \mathbf{R}}=\underset{t=1}{\overset{{{N}_{i}}}{\sum}}\,\frac{\partial {{C}_{t}}}{\partial \mathbf{R}}\end{eqnarray}$$

C_total and $\frac{\partial {{C}_{\text{total}}}}{\partial \mathbf{R}}$ can then be used to find the optimal values of R using the same gradient based optimization methods as used for standard registration between two images.

With this new framework it is possible to fit the model directly to the partial image data, removing the need for full images, by modelling the image acquisition process:

$$\begin{eqnarray}&&{{\mathbf{P}}_{t}}={{A}_{t}}\left({{\mathbf{I}}_{t}}\right)+\boldsymbol{\varepsilon }_{{t}}\end{eqnarray}$$

where P_t is the measured partial image data at time t, I_t the (unknown) full image, and ε_t the imaging noise. A_t is the function which simulates the image acquisition, e.g. for projection data A_t would be the forward-projection operator and for slice data A_t would be the slice selection profile.

When fitting the correspondence model the partial image data is simulated from the transformed reference-state image, and the similarity term, Sim, measures the (dis)similarity between the measured and simulated partial image data, i.e.

$$\begin{eqnarray}&&{{C}_{t}}\left({{\mathbf{P}}_{t}},{{\mathbf{I}}_{0}},{{\mathbf{M}}_{t}}\right)=Sim\left({{\mathbf{P}}_{t}},{{A}_{t}}\left(T\left({{\mathbf{I}}_{0}},{{\mathbf{M}}_{t}}\right)\right)\right)+Con\left({{\mathbf{M}}_{t}}\right)\end{eqnarray}$$

However, Sim, is now defined in the space of the partial image data but the transformation parameterized by M_t is defined in the space of the full images. Therefore, when calculating $\frac{\partial {{C}_{t}}}{\partial {{\mathbf{M}}_{t}}}$ it is necessary to transform the gradient of Sim from partial image space into full image space. This is done using the adjoint of the image acquisition function, $A_{t}^{\ast}$, e.g. the adjoint of the forward projection operator is the back projection operator:

$$\begin{eqnarray}&&\frac{\partial {{C}_{t}}}{\partial {{\mathbf{M}}_{t}}}=\frac{\partial Sim}{\partial {{\mathbf{M}}_{t}}}+\frac{\partial Con}{\partial {{\mathbf{M}}_{t}}}\end{eqnarray}$$

$$\begin{eqnarray}&&\frac{\partial Sim}{\partial {{\mathbf{M}}_{\boldsymbol{t}}}}=\frac{\partial {{\mathbf{I}}_{{{T}_{t}}}}}{\partial {{\mathbf{M}}_{\boldsymbol{t}}}}\frac{\partial Sim}{\partial {{\mathbf{I}}_{{{T}_{t}}}}}\end{eqnarray}$$

$$\begin{eqnarray}&&\frac{\partial Sim}{\partial {{\mathbf{I}}_{{{T}_{t}}}}}=A_{t}^{\ast}\left(\frac{\partial Sim}{\partial {{\mathbf{P}}_{{{A}_{t}}}}}\right)\end{eqnarray}$$

where ${{\mathbf{I}}_{{{T}_{t}}}}=T\left({{\mathbf{I}}_{0}},{{\mathbf{M}}_{t}}\right)$, and ${{\mathbf{P}}_{{{A}_{t}}}}={{A}_{t}}\left({{\mathbf{I}}_{{{T}_{t}}}}\right)$. $\frac{\partial Sim}{\partial {{\mathbf{P}}_{{{A}_{t}}}}}$ is the gradient of the similarity measure with respect to the simulated partial image data, $\frac{\partial Sim}{\partial {{\mathbf{I}}_{{{T}_{t}}}}}$ is the gradient of the similarity measure with respect to the transformed reference-state image, and $\frac{\partial {{\mathbf{I}}_{{{T}_{t}}}}}{\partial {{\mathbf{M}}_{t}}}$ is the gradient of the transformed reference-state image with respect to the motion parameters, e.g. for B-spline transformations this is the spatial gradient of the reference-state image transformed by M_t and convolved with the B-spline kernel.

The framework as described above assumes a reference-state image, I₀, is available. However, if the motion is known MCIR can be used to reconstruct I₀ from the partial image data (Batchelor et al 2005, Rit et al 2009, Van Reeth et al 2012), so a separate reference-state image is not required. MCIR can be incorporated into this framework using an iterative approach. Firstly, a MCIR is performed from the partial image data using an initial estimate of no motion (or using a previous model if one is available, e.g. from an earlier scan). The result will contain blurring and/or other motion artefacts, but can be used as an initial estimate of I₀ to fit the model parameters, R, as described above. The fitted model is then used to perform another MCIR, updating the estimate of I₀. The process continues to iterate between fitting R and reconstructing I₀ until there is no more improvement or the maximum number of iterations is reached.

The framework has been implemented using the open source NiftyReg⁷ deformable image registration software (Modat et al 2010). This uses a B-spline transformation model and conjugate gradient optimisation in a multi-resolution approach. As all the experiments use intra-modality data, sum of squared differences (SSD) was used as the similarity measure. The use of a constraint term was found to be unnecessary.

Any correspondence model can be used which can be expressed as a linear combination of terms pre-calculated from the surrogate signal(s), ${{\varphi}_{n}}\left({{\mathbf{S}}_{t}}\right)$, i.e.

$$\begin{eqnarray}&&{{\mathbf{M}}_{t}}=\underset{n=1}{\overset{{{N}_{r}}}{\sum}}\,{{\mathbf{R}}_{n}}{{\varphi}_{n}}\left({{\mathbf{S}}_{t}}\right)\end{eqnarray}$$

then

$$\begin{eqnarray}&&\frac{\partial {{\mathbf{M}}_{t}}}{\partial {{\mathbf{R}}_{n}}}={{\varphi}_{n}}\left({{\mathbf{S}}_{t}}\right)\end{eqnarray}$$

This includes linear, polynomial, and B-spline models. e.g. for a linear model with two surrogate signals, ${{\mathbf{S}}_{t}}=\left[{{s}_{t,1}},{{s}_{t,2}}\right]$:

$$\begin{eqnarray}&&{{\varphi}_{1}}\left({{\mathbf{S}}_{t}}\right)={{s}_{t,1}},\quad {{\varphi}_{2}}\left({{\mathbf{S}}_{t}}\right)={{s}_{t,2}},\quad {{\varphi}_{3}}\left({{\mathbf{S}}_{t}}\right)=1\end{eqnarray}$$

for a 2nd order polynomial model with one surrogate signal, ${{\mathbf{S}}_{t}}={{s}_{t}}$:

$$\begin{eqnarray}&&{{\varphi}_{1}}\left({{\mathbf{S}}_{t}}\right)=s_{t}^{2},\quad {{\varphi}_{2}}\left({{\mathbf{S}}_{t}}\right)={{s}_{t}},\quad {{\varphi}_{3}}\left({{\mathbf{S}}_{t}}\right)=1\end{eqnarray}$$

and for a periodic B-spline model with 4 control points (N_r  =  4) and respiratory phase as the surrogate signal, ${{\mathbf{S}}_{t}}={{\vartheta}_{t}}$ (where ${{\vartheta}_{t}}$ is the respiratory phase and has a value between 0 and 1):

$$\begin{eqnarray}&&{{\varphi}_{n}}\left({{\mathbf{S}}_{t}}\right)={{B}_{i}}\left(\,j\right)\end{eqnarray}$$

Where $i=\left(\lfloor 4{{\vartheta}_{t}}\rfloor +n-2\right)\,\bmod \,4$, $j=4{{\vartheta}_{t}}-\,\lfloor 4{{\vartheta}_{t}}\rfloor$, and ${{B}_{i}}$ is the ith cubic B-spline basis function:

$$\begin{eqnarray}&&{{B}_{1}}\left(\,j\right)={{(1-j)}^{3}}/6\end{eqnarray}$$

$$\begin{eqnarray}&&{{B}_{2}}\left(\,j\right)=\left(3{{j}^{3}}-6{{j}^{2}}+4\right)/6\end{eqnarray}$$

$$\begin{eqnarray}&&{{B}_{3}}\left(\,j\right)=\left(-3{{j}^{3}}+3{{j}^{2}}+3j+1\right)/6\end{eqnarray}$$

$$\begin{eqnarray}&&{{B}_{4}}\left(\,j\right)={{j}^{3}}/6\end{eqnarray}$$

The dynamic images can be full images, 'slab' images (consisting of a few slices), or individual slices. If the dynamic images are not full images the image acquisition function, A_t, samples the deformed image, ${{\mathbf{S}}_{{{T}_{t}}}}$, at the slice(s) corresponding to the dynamic image, P_t. The dynamic images can have a different resolution to the reference-state image, I₀, enabling the use of a high resolution reference-state image with low resolution dynamic images, and permitting the use of super-resolution MCIR methods. If the dynamic images are lower resolution than I₀ the function A_t first convolves the deformed image, $\boldsymbol{I}_{{{{T}_{t}}}}$ (which has the same resolution as I₀), with an appropriate Gaussian kernel to account for the different image resolutions (Cardoso et al 2015), before sampling the slices corresponding to P_t. Likewise, the adjoint function, $A_{t}^{\ast}$, first resamples $\frac{\partial Sim}{\partial {{\mathbf{P}}_{{{A}_{t}}}}}$ from the space of ${{\mathbf{P}}_{t}}$ into the space of the full deformed image, ${{\mathbf{I}}_{{{T}_{t}}}}$. If the dynamic images are lower resolution, the result is then convolved with the same Gaussian kernel as above.

Two MCIR methods have been implemented: a simple averaging of the deformed dynamic images, and a super-resolution reconstruction using the iterative back-projection method (Irani and Peleg 1991). When performing a MCIR the the dynamic image data (after applying $A_{t}^{\ast}$, i.e. $A_{t}^{\ast}\left({{\mathbf{P}}_{t}}\right)$) must be deformed into the space of I₀. This could be done by estimating the inverse of the transformations parameterized by M_t, but this can be computationally expensive. Instead the adjoint function T ^* is used, i.e. 'push-interpolation' is used instead of the usual 'pull-interpolation'. This can potentially lead to 'holes' (voxels with no intensity information) in the deformed images, but as long as the holes are not at the same location for all dynamic images they will not be present in the final MCIR. If holes are present in the final MCIR, $A_{t}^{\ast}\left({{\mathbf{P}}_{t}}\right)$ can be resampled at a higher resolution before applying T^*. However, sampling $A_{t}^{\ast}\left({{\mathbf{P}}_{t}}\right)$ at the same resolution as I₀ and using linear interpolation resulted in no holes in the MCIRs for all the experiments below.

Four real data sets have been used: full images (4DCT), slab images (Cine CT), thin slices (helical CT), and thick slices (multi-slice MR).

The full images dataset is one of the freely available DIR-Lab 4DCT datasets⁸ (Castillo et al 2009). This consists of 10 CT volumes with a resolution of 1.1  ×  1.1  ×  2.5 mm³, representing different phases of the respiratory cycle.

The slab image dataset uses Cine CT volumes (unsorted 4DCT data) acquired at University College Hospital, London, as the dynamic images. Each Cine CT volume consists of 8  ×  2.5 mm slices (0.98 mm  ×  0.98 mm in-slice resolution), i.e. each volume only covers 2 cm, but they are acquired from 17 contiguous couch-positions such that the entire lungs are covered. Each couch-position is imaged 11 times, sampling just over one breath cycle, giving a total of 187 dynamic images. The surrogate signal was acquired using the Varian RPM system.

The thin slice dataset consists of individual CT slices from fast helical scans acquired at UCLA (Thomas et al 2014). Each helical scan consists of 301 slices and takes ~5 s to acquire. The images were resampled with a resolution was 1  ×  1  ×  1 mm³. 25 helical scans were acquired in alternating directions while the patient was freely breathing, although only 10 scans were used here. Therefore, each individual slice is imaged 10 times, with the images corresponding to arbitrary points in different breath cycles, such that both the intra- and inter-cycle variations are sampled for each slice. The surrogate signal was acquired using a respiratory belt.

The thick slice dataset consists of multi-slice MRI data acquired at the Institute of Cancer Research, London. Sagittal multi-slice MR data were acquired from a volunteer during free-breathing, from 30 adjacent slices encompassing both lungs. Each slice was imaged 10 times. The in-slice resolution was 1.77 mm  ×  1.77 mm and the slices thickness was 10 mm. The acquisition time for each slice was approximately 0.2 s, so the full acquisition took approximately 1 min. This acquisition was then repeated 4 more times, with a 2 mm left-right shift between each acquisition, giving a total of 150 overlapping slices. The surrogate signal was acquired from an MR-compatible ABC device working passively as a spirometer (Kaza et al 2015).

The details of the experiments performed on the different datasets are given in table 1. For the full images dataset a periodic B-spline correspondence model was used with respiratory phase as the surrogate signal. This can model intra-cycle variation (hysteresis), but not inter-cycle variation. For the other datasets a linear correspondence model was used, relating the motion to both the value and temporal derivative of the surrogate signal, i.e. ${{\mathbf{S}}_{t}}=\left[{{s}_{t}},{{\overset{\centerdot}{{s}}\,}_{t}}\right]$. This correspondence model has been popular in the literature due to its simplicity and ability to model both intra- and inter-cycle variations. The control point grid (CPG) spacing of the B-spline transformation was empirically set for each experiment, such that reasonable results were obtained. In all cases, the same spacing was used in all directions.

For experiment 1 on the full images and slab images datasets the end-exhalation 4DCT volumes was used as I₀, as this is considered the most reproducible position and so the images usually contain the fewest artefacts. For experiment 1 on the thin slice dataset one of the helical scans was used as I₀. As this was acquired while the patient was freely breathing the anatomy appears distorted, but due to the fast scan time there were no obvious missing or repeated structures. For the experiments that used the simple averaging MCIR method I₀ was reconstructed with the same resolution as the dynamic images. For experiment 2 on the thick slice dataset the super-resolution MCIR method was used to reconstruct I₀ with a voxel size of 2  ×  1.77  ×  1.77 mm³.

Note, as the surrogate signals have been normalized so that their mean values equal 0 the reference-state images reconstructed by MCIR will represent the average position of the anatomy during the acquisition. For the experiments that did not perform MCIR, I₀ does not represent the average position of the anatomy. Therefore, and extra constant offset model parameter was included in the linear models to account for the difference between the I₀ and the average position of the anatomy (this is not necessary for the B-spline model).

One problem encountered with several of the datasets was sliding motion between the lungs and the chest wall, which cannot be correctly modelled by a standard B-spline transformation. To alleviate this problem, the lungs were automatically segmented in some datasets using thresholding and morphological operations, as indicated in table 1. Voxels outside the lungs were set to the same intensity as soft tissue. For the slab images automatic lung segmentation was challenging in the dynamic images. Therefore, for experiment 1 the lungs were only segmented in I₀, and for experiment 2, which used MCIR, they were not segmented at all. For the thin slice data good results were obtained with unsegmented images, so segmentation was not used.

For the full images manually located landmark points are available which can be used to quantitatively assess the results. These consist of 300 landmarks well distributed over both lungs, which have been located in the end-inhalation and end-exhalation images. A subset of 75 of the landmarks have also been located on all of the exhalation images. Both sets of landmarks were used to assess the fitted motion model using the 'snap-to-voxel' approach described in Castillo et al (2009). For comparison the landmark errors have also been calculated using an estimate of no motion, i.e. measuring the magnitude of the motion in the original images.

For the other datasets it is difficult to accurately locate landmarks in the images, therefore only qualitative assessments of their results have been possible.

As it is difficult to estimate the ground truth motion for most of the real datasets, additional experiments were performed using datasets generated from a computer phantom. For these experiments the ground truth motion and reference-state image, I_true, are known, so the results of the motion model fitting and the MCIR can be quantitatively assessed. A simple 2D 'lung-like' phantom was used for all the experiments (figure 2(a)). The phantom was animated using displacement fields generated from a linear correspondence model relating the motion to a surrogate signal and its temporal gradient. For the slab images, thin slices, and thick slices, the surrogate signals from the real datasets were used to animate the phantom. For the full images a real surrogate signal is not available, so an idealized surrogate signal from a cos⁴ function was used. For all of the phantom image datasets Gaussian noise with a standard deviation of 3% of the maximum image intensity (1500) was added to the simulated images. When animating the phantom, the displacement field itself was used as the motion parameters, M_t, whereas the B-spline control point displacements are used as M_t when fitting the model. This means there are two motion parameters at each pixel, the displacement in the x and y directions, and therefore four correspondence model parameters for each pixel, r_x,1, r_x,2, r_y,1, and r_y,2. The values of these model parameters were manually defined such that they were smoothly varying over the phantom (figure 2(b)), and produced plausible looking motion that included deforming anatomy and intra- and inter-cycle variations.

Zoom In Zoom Out Reset image size

The phantom datasets and the experiments performed with them mimicked the real datasets and experiments as closely as possible. The same number of images/slabs/slices was used as for the real data, but as the phantom was only 2D, images were 2D, and slices were 1D. Different sized phantoms were required for the different datasets: 128  ×  128 pixels for the full images, 136  ×  136 pixels for the slab images, 301  ×  301 pixels for the thin slices, and 150  ×  150 pixels for the thick slices. However, the phantom was resized so that its proportions were always the same. As the phantom used for the thin slices had approximately twice as many pixels as for the other datasets, the magnitude of motion was twice as large (in pixels) for this dataset. The thick slices were simulated by convolving the deformed images with a Gaussian kernel representing a slice thickness of 5 pixels. The surrogate signals from the real datasets were also used for the phantom experiments, and the same correspondence models and MCIR methods were applied as for the real data. The B-spline CPG spacing was empirically set to 10 pixels for the thin slice dataset, and to 5 pixels for all the others. As the phantom datasets did not contain sliding motion there was no need for segmentation.

The results of all the phantom experiments were assessed by calculating the Displacement Field Error, DFE, defined as the 2D Euclidean distance at each pixel between the true displacement field and the displacement field resulting from the fitted motion model. The DFE was calculated at every pixel inside the deformed phantom for every time-point. Pixels outside the phantom were ignored as they do not contain any image data to guide the model fitting. The DFE was also calculated for each experiment using an estimate of no motion, to quantify the amount of motion present in the original data. The experiments that performed MCIRs were also assessed by comparing the reconstructed I₀ to I_true. This was accomplished by calculating the absolute difference in the image intensities at every pixel and the correlation coefficient between the images. Similar calculations were performed for an image reconstructed using the same method (averaging or super-resolution) but assuming no motion.

Videos showing the results of the experiments on the real datasets are available as supplementary files. These show the original dynamic images, I_i or P_i, the simulated dynamic images, ${{\mathbf{I}}_{{{T}_{i}}}}$ or ${{\mathbf{P}}_{{{A}_{t}}}}$, and a colour overlay of these. In addition, for the datasets using partial image data (i.e. the slab and slice datasets), the videos also show one or more views through the full deformed reference-state image, ${{\mathbf{I}}_{{{T}_{i}}}}$.

4.1.1. Full images.

4.1.2. Slab images.

There was a fairly good match between P_i and ${{\mathbf{P}}_{{{A}_{t}}}}$ for both experiments, as can be seen in the supplementary movies (real_data_slab_images_exp1.mp4 and real_data_slab_images_exp2.mp4). As P_i only sample one respiratory cycle at each couch-position it is difficult to assess how accurately the models estimate the inter-cycle variations, but the images of ${{\mathbf{I}}_{{{T}_{i}}}}$ show plausible motion over the entire field of view, exhibiting both intra- and inter-cycle variations. As with the 4DCT data, there are some noticeable mismatches at the back of the lungs due to the sliding motion, and some smaller mismatches in other regions, but again these are partly due to blurred and repeated structures in some of the Cine CT volumes.

Figure 3 compares the end exhalation 4DCT volume with the MCIR produced in experiment 2, and a similar image formed without applying any motion correction (i.e. the result of averaging all of the Cine CT volumes together). It is evident that blurred structures are greatly reduced in the MCIR compared to the result with no motion correction, indicating that most of the motion has been well recovered. Some structures still appear slightly blurred compared to the 4DCT volume, particularly near the back of the lungs where sliding motion occurs. This blurring is partly caused by errors in the estimated motion, but also due to the structures appearing blurred in some of the Cine CT volumes.

Zoom In Zoom Out Reset image size

4.1.3. Thin slices.

For both experiments there is a very good match between P_i and ${{\mathbf{P}}_{{{A}_{t}}}}$, as can be seen in the supplementary movies (real_data_thin_slices_exp1.mp4 and real_data_thin_slices_exp2.mp4). The images of ${{\mathbf{I}}_{{{T}_{i}}}}$ show plausible motion, including intra- and inter-cycle variations, over the entire field of view, and even seem to reproduce most of the sliding motion seen at the back of the lung. Since P_i sample both the intra- and inter-cycle variations at each location, this implies the true variations have been modelled well. Figure 4 shows the helical scan used as I₀ for experiment 1, the MCIR produced by experiment 2, and a similar image formed without applying any motion correction. The blurring has been almost completely removed in the MCIR, with the majority of structures appearing as sharp as in the helical scan, indicating the motion has been modelled very well. In addition, the signal-to-noise ratio was noticeably improved in the MCIR due to combining the image data from all 10 slices at each location.

Zoom In Zoom Out Reset image size

4.1.4. Thick slices.

${{\mathbf{P}}_{{{A}_{t}}}}$ match P_i well for both experiments, as can be seen in the supplementary movies (real_data_thick_slices_exp1.mp4 and real_data_thick_slices_exp2.mp4). Similar to the thin slices dataset, P_i sample both the intra- and inter-cycle variations at each location, which indicates these have been modelled well. The images of ${{\mathbf{I}}_{{{T}_{i}}}}$ show plausible motion over the entire field of view, although the motion at the edge of the lungs looks questionable, particularly for experiment 2. This is due to the sliding motion between the lungs and the ribs, and the ribs sometimes being included in the segmented lungs as they have similar intensity in the original MR images.

Figure 5 shows an example P_i (after lung segmentation), together with the corresponding ${{\mathbf{P}}_{{{A}_{t}}}}$ from both experiments. The signal-to-noise ratio has greatly improved in ${{\mathbf{P}}_{{{A}_{t}}}}$ from experiment 1, due to averaging the dynamic images, and ${{\mathbf{P}}_{{{A}_{t}}}}$ from experiment 2 shows even finer structure inside the lung due to the super-resolution reconstruction and an increased number of dynamic images. Figure 6 shows coronal views through the MCIRs produced by both experiments, together with similar images formed without applying any motion correction. The blurring has been removed in the MCIRs, indicating that motion has been recovered well. Much more detail can be seen in the MCIR from experiment 2, indicating that the super-resolution method has been applied successfully and that the motion has been modelled well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Videos showing the results of the experiments on the phantom datasets are available as supplementary files. Each of these show the original dynamic images, I_i or P_i, the simulated dynamic images, ${{\mathbf{I}}_{{{T}_{i}}}}$ or ${{\mathbf{P}}_{{{A}_{t}}}}$, and a colour overlay of these. Also, for the phantom datasets using partial image data the true deformed reference-state images, I_i, are shown together with the reference-state image deformed using the fitted model, ${{\mathbf{I}}_{{{T}_{i}}}}$, and a colour overlay of these. The videos show that there were good matches between P_i and ${{\mathbf{P}}_{{{A}_{t}}}}$ and between I_i and ${{\mathbf{I}}_{{{T}_{i}}}}$ for all phantom datasets experiments. In addition, for the experiments that performed MCIRs, I₀ closely resemble I_true, although for experiment 1 on the thick slice data I₀ has a noticeably lower through-slice resolution, as expected.

The qualitative results shown in the videos are supported by the quantitative results in tables 3 and 4. Table 3 gives the Displacement Field Error (DFE) for each experiment. When using the fitted model the DFE is greatly reduced for all experiments, in most cases to below 1 pixel, indicating the motion has been recovered very well. The DFE is larger for experiment 1 on the thin slice data, but the magnitude of the motion was twice as large for this dataset (as it represents higher resolution data), and the relative reduction in DFE when using the model is comparable to the other datasets. Interestingly, the use of MCIR (i.e. experiment 2) leads to a considerably lower DFE for the thin slice dataset but not for the slab images dataset. For the thick slice data, experiment 2, which uses a super-resolution MCIR, yields a much lower DFE than experiment 1, which uses the averaging MCIR. This result is expected, as the super-resolution MCIR has a higher resolution and hence can recover the motion more accurately.

Table 4 displays the results of comparing I₀ to I_true, and figure 7 shows images of I₀ for each experiment performing MCIR. These results indicate that the motion has been recovered well by the fitted motion models and that I₀ closely resembles I_true for all experiments. As expected, the results for the thick slice dataset are not as good as for the thin slice and slab image datasets, because the dynamic images have a lower resolution than I_true. However, the use of the super-resolution method improves the results, and produces an I₀ more similar to I_true.

Zoom In Zoom Out Reset image size