MoVIT: a tomographic reconstruction framework for 4D-CT

Vincent Van Nieuwenhove; Jan De Beenhouwer; Jelle Vlassenbroeck; Mark Brennan; Jan Sijbers

doi:10.1364/OE.25.019236

1. Introduction

4D (3 spatial dimensions + time) computed tomography (CT) has the ability to image dynamic processes in a non-destructive manner. It has great potential for validating models of dynamic processes, such as the compression of polyurethane (PU) foam, biomechanical properties of organs, baking or leavening of pastry products [1–3].

In 4D-CT, several CT datasets (time frames) are typically acquired sequentially [1]. Each time frame consists of a set of projections acquired over an angular range of 180 or 360 degrees. Afterwards, these time frames are independently reconstructed. However, due to the long acquisition time of conventional CT, two problems arise if a fast dynamic process is imaged. Firstly, the long acquisition time of a single time frame strongly limits the temporal resolution. Secondly, due the object deformation during the acquisition of a single time frame, the projections of this time frame are not consistent with each other. This results in blurry reconstructed images due to deformation artefacts.

A straightforward method to avoid both problems is shortening the acquisition time of a single time frame. This can be achieved by reducing the integration time of a single projection or lowering the number of projections per time frame. A shorter integration time, however, reduces the signal-to-noise ratio (SNR) of the projection data, which in turn leads to a lower SNR in the reconstructed images. Lowering the number of projections, on the other hand, results in streaks in the reconstructed images. As a result, the conventional workflow leads to a trade-off between the temporal resolution/deformation artefacts and low SNR/streaking artefacts in the reconstructed images.

Fortunately, the trade-off between the temporal resolution and SNR can be improved by exploiting data redundancy present in 4D-CT datasets. Since in every time frame the same, though slightly changed, object is scanned, it is beneficial to include information about other time frames into the reconstruction process [4]. In this paper we will focus on non-affine deformation correction, in contrast with several proposed algorithms [5–8]. Previously proposed reconstruction methods exploiting a non-affine motion model can be split up into different classes. A first class assumes structural changes to be local. Such algorithms either rely on an a priori high quality reconstruction or assume that the location of structural changes is known [4,9,10]. A second class regularises both the spatial and temporal domain of the reconstructions with, for example, Markov random fields [11,12]. While these algorithms do not depend on a particular deformation model they often suffer from long computation times and a large number of tunable parameters.

The last class of methods assumes the deformation to be described as a diffeomorphic deformation vector field (DVF) that is included into the reconstruction algorithm. This vector field describes the displacement of every voxel between two time frames. Since the deformation is a priori unknown, it has to be estimated. Some reconstruction methods include volume registration to improve reconstructions with deformation artefacts due to respiratory and cardiac motion [13–16]. These methods exploit the periodicity of the motion [13], assume that a prior deformation free scan is available [15], assume multiple quiescent motion time frames [14], or have a high number of tunable parameters [16]. These properties prevent the use of these algorithms in the context of micro-CT where a wide range of samples and dynamics has to be studied.

In this paper, we propose the MotionVector-based Iterative Technique (MoVIT). It incorporates the deformation fields directly into the reconstruction process and, as such, exploits information from projections of other time frames into the reconstruction of the image at a particular time frame. In order to estimate the deformation field map, non-rigid volume registrations are performed on conventionally reconstructed volumes of the different time frames.

The paper is structured as follows. In section 2, the proposed reconstruction algorithm and deformation estimation strategy, MoVIT, is explained. Next, in section 3, MoVIT is compared to common reconstruction algorithms in a numerical way as well as on an experimental dataset, the results of which are discussed in section 4.

2. Method

In subsection 2.1, an overview of the stationary tomographic model will be given and generalised (in subsection 2.2) to the dynamic case. Next, in subsection 2.3, the MoVIT reconstruction algorithm is described, which incorporates a deformation field into an algebraic reconstruction method (ARM) that is assumed to be known a priori. This restriction is then abrogated in subsection 2.4, where an estimation method for the deformation vector field is introduced.

2.1. Stationary algebraic tomography

Let x = (x_i) ∈ ℝ^N be a vector of unknown volume elements representing the scanned object. The log-corrected projection values for all projection angles θ = (θ_i) ∈ [0, 2π[^M are denoted by p = (p_i) ∈ ℝ^M. The projection process can be simulated as q = Ax, where A = (a_ij) ∈ ℝ^M×N is a matrix of which the entries a_ij represent the contribution of pixel value x_j to projection value q_i (see Fig. 1).

Fig. 1 Illustration of the projection process. The contribution a_ij of pixel x_j to the projection value q_i is represented as the ray-intersection length of projection line i with pixel j.

Download Full Size | PDF

The system of linear equations Ax = p cannot be solved directly since A is generally not invertible. A closed form expression for the (regularized) least-squares solution can be derived. However, due to the size of the problem, the direct calculation of this solution is infeasible on modern computers. Therefore, algebraic reconstruction methods such as ART, SART or SIRT, start with an initial guess x = x⁰ and iteratively compute new estimates x^k (k = 1, 2, . . .). This is repeated until an approximate solution of Ax = p is found. For example, the simultaneous iterative reconstruction technique (SIRT) algorithm solves the weighted least-squares optimization problem argmin_x ‖Ax − p‖_R, where ‖Ax − p‖_R = (Ax − p)^T R(Ax − p) and R = (r_kl) ∈ ℝ^M×M is a diagonal matrix with $r_{k k} = {(\sum_{l} a_{k l})}^{- 1}$ [17]. The following iterative formula is known to converge to the weighted least-squares minimum:

x^{k} = x^{k - 1} + {CA}^{T} R (p - {Ax}^{k - 1})

where C = (c_kl) ∈ ℝ^N×N is a diagonal matrix with

c_{l l} = {(\sum_{k} a_{k l})}^{- 1}

.

2.2. Dynamic tomographic model

The conventional algebraic tomography model described in section 2.1 assumes the scanned object to remain stationary throughout the acquisition process. This assumption is no longer valid in dynamic CT. Therefore, the standard model has to be generalized to deal with dynamic objects.

A dynamic object can be represented as a time series of images x_r ∈ ℝ^N, where r ∈ {1, . . . , R} is the time index, with R the total number of time frames. During the acquisition, the gantry rotates multiple times around the object while acquiring projections. The projections of subscan r are represented by p_r ∈ ℝ^M_r. The sparse matrix A_r ∈ ℝ^{M_r ×N} is the corresponding forward projection matrix. If the object is assumed stationary during each time frame, the acquisition of the dynamic process can be modelled as follows:

[\begin{matrix} A_{1} & 0 & \dots & 0 \\ 0 & A_{2} & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & A_{R} \end{matrix}] = [\begin{matrix} x_{1} \\ x_{2} \\ ⋮ \\ x_{R} \end{matrix}] = \tilde{A} \tilde{x} = \tilde{p},

where Ã represents the block diagonal matrix consisting of blocks A₁, . . . , A_R,

\tilde{x} = {(x_{1}^{T}, \dots, x_{R}^{T})}^{T} \in ℝ^{R N}

and

\tilde{p} = {(p_{1}^{T}, \dots, p_{R}^{T})}^{T} \in ℝ^{\sum_{r = 1}^{R} M_{r}}

[4, 11]. Notice that the dynamic computed tomography model does not take deformation during a single time frame into account. As a result, deformation artefacts will occur if excessive motion takes place during a single time frame (more than a single voxel). To remedy this, the acquisition time of a single time frame can be reduced.

The dynamic model (2), however, does not include the dependencies of the object between different time frames. The temporal dependencies between the subsequent images x_r can be modelled as follows. Each image x_r is regarded as a discretisation of f_r (y) (y ∈ ℝ³). A deformation between two time points can be described as; f_r (y) = f_r′ (y+v_r′r (y)) with v_r′r (y) the DVF from time point r′ to r. The discretized version of the previous equation comes down to interpolation:

x_{r} = τ_{r^{'} r} x_{r^{'}}

with the operator τ_r′r : ℝ^N → ℝ^N transforming the image of the object at time frame r′ to its state at time frame r. For many interpolation methods (e.g., nearest neighbour, linear or cubic interpolation), the operator τ_r′r is actually a linear operator and thus τ_r′r ∈ ℝ^N×N.

The transformations τ_·· can be included in the tomographic model:

\tilde{A} \tilde{x} = \tilde{A} [\begin{matrix} x_{1} \\ ⋮ \\ x_{R} \end{matrix}] = \tilde{A} [\begin{matrix} τ_{r 1} \\ ⋮ \\ τ_{r R} \end{matrix}] x_{r}

and thus,

𝒜_{r} x_{r} = \tilde{p}

with

𝒜_{r} = [\begin{matrix} A_{1} τ_{r 1} \\ ⋮ \\ A_{R} τ_{r R} \end{matrix}],

the operator transforming the object at time frame r to the projections p̃, which are acquired over the whole experiment (or a subset of the time frames). Notice that in (5) the time dependency is modelled by the operator 𝒜_r, in contrast with (2). As a result, the number of unknowns is reduced from RN to N with the same number of equations.

Note that modelling the temporal dependencies with DVFs allows for an accurate description if the deformation is a diffeomorphism. In other words, the DVFs are differentiable and have a differentiable inverse, such as elastic deformations. Unfortunately, fluid flow and structural changes such as the formation of cracks cannot be accurately described with these models. For fluid flow specialised reconstruction algorithms do exist [5,11]. In future research these algorithms could be combined with the proposed MoVIT framework to take both kinds of deformations into account.

2.3. Motion vector-based iterative technique

Assuming that the deformation of the object during the acquisition is known, an approximation of the solution of (5) can be found with the MoVIT algorithm. The basic steps of MoVIT are shown in Algorithm 1. In iteration k, the current reconstruction of time frame r, $x_{r}^{k}$ , is transformed to each of the other time frames r′ ∈ 𝒩_r with the deformation operator τ_rr′ and 𝒩_r the set of neighbouring time frames of time frame r. With this transformed reconstruction and the projections of that particular time frame, an ARM (e.g., ART, SART, SIRT, ...) update q is calculated. These ARM updates are then transformed back to the time frame r using $τ_{r r^{'}}^{- 1}$ with weight w_rr′. The weights w_rr′ are normalized such that $\sum_{r^{'} \in 𝒩_{r}} w_{r r^{'}} = 1$ . The proposed calculation of the weights will be explained in Section 2.4. Afterwards, the transformed updates added to the MoVIT update u. This MoVIT update is then added to the current reconstruction, resulting in $x_{r}^{k + 1}$ . This scheme is repeated until the stop criterion is met or a maximum number of MoVIT iteration is reached.

Algorithm 1. Basic steps of the MoVIT algorithm

View Table | View all tables in this article

The MoVIT framework can be implemented with any algebraic reconstruction method such as ART, SART, or SIRT. In the rest of this work, SIRT will be the method underlying MoVIT [17]. To reconstruct the object at time frame r, the framework results in following algorithm:

x_{r}^{k + 1} = x_{r}^{k} + \sum_{r^{'} \in 𝒩_{r}} w_{r r^{'}} τ_{r r^{'}}^{- 1} C_{r^{'}} A_{r^{'}}^{T} R_{r^{'}} (p_{r^{'}} - A_{r^{'}} τ_{r r^{'}} x_{r}^{k})

where C_r = (c_r,kl) ∈ ℝ^N×N is a diagonal matrix with

c_{r, l l} = {(\sum_{k} a_{k l})}^{- 1}

and R_r = (r_r,kl) ∈ ℝ^{M_r ×M_r} is a diagonal matrix with

r_{r, k k} = {(\sum_{l} a_{k l})}^{- 1}

. An overview of a single MoVIT iteration is given in Fig. 2.

Fig. 2 Schematic representation of a single iteration of the MoVIT algorithm implemented with SIRT updates.

Download Full Size | PDF

2.4. Deformation estimation

The deformation operators τ_rr′ and their inverse (see Algorithm 1) are unknown and have to be estimated. To this end, each time frame is first conventionally reconstructed, for example with SIRT, using only the projections corresponding to that time frame. Afterwards, these conventional reconstructions (c₁, c₂, . . . , c_R) are pairwise registered with each other, resulting in DVFs. In this work, the registration needed to estimate the parameters μ_rr′ of a b-spline deformation model was performed with Elastix [18]. The metric of the image registration algorithm is the mean squared difference (MSD) and thus the registration algorithm solves:

μ_{r r^{'}} = {argmin}_{μ} (\frac{1}{N} \sum_{i = 1}^{N} {(c_{r^{'}, i} - {(τ (μ) c_{r})}_{i})}^{2})

where c_r′,i is the i^th element of c_r′ and (τ(μ)c_r)_i is the i^th element of τ(μ)c_r. The transform τ_r′r is in theory also the inverse of τ_rr′. However, the exact solution of non-rigid image registration is often non-unique or nonexistent. Therefore, it is crucial to determine the direct inverse deformation field of the obtained DVF. This inverse DVF can be calculated with the method described in Chen et al. [19]. Evidently, it is unnecessary to determine the transform τ_rr and its inverse as they are the identity transform. Note that other image registration algorithms, such as optical flow algorithms and digital volume correlation methods, are compatible as well with the MoVIT frame work [20]. The weights w_rr′ reflect the accuracy of the corresponding DVFs and are calculated as follows:

w_{r r^{'}} = \frac{exp (- {(k_{r r^{'}} / b)}^{2})}{\sum_{l \in 𝒩_{r}} exp (- {(k_{r l} / b)}^{2})}

where

k_{r r^{'}} = N^{- 1} \sum_{i} {(c_{r, i} - {(τ (μ_{r r^{'}}) c_{r}^{'})}_{i})}^{2}

corresponds to the MSD between the reconstructed time frame and the transformed reconstruction of the time frame r′. The parameter b regulates the magnitude of the weights. If the MSD of a time frame equals b, the weight of that time frame will be exp(1) times smaller than the weight of the reconstructed time frame. The proposed method was implemented in Matlab, and the forward and back projections were performed with the ASTRA toolbox [21–23]. An overview of the complete framework is shown in Fig. 3.

Fig. 3 Schematic representation of the full MoVIT reconstruction pipeline.

Download Full Size | PDF

3. Experiments and results

In subsection 3.1, the reconstruction methods, employed in the experiments, are explained. The MoVIT method is compared to conventional reconstruction methods in numerical simulations (subsection 3.3) and on an experimental dataset of polyurethane foam under compression (subsection 3.4). The results of these experiments are discussed in section 4.

3.1. Reconstruction methods

The proposed MoVIT algorithm was compared against two conventional reconstruction algorithms which are independently applied on each time frame. The first conventional algorithm is the Feldkamp-David-Kress (FDK) algorithm, resulting in the reconstructions f₁, f₂, . . . , f_R, where f_r ∈ ℝ^N [24]. The second algorithm is the SIRT algorithm (see Section 2.1) which results in the reconstructions c₁, c₂, . . . , c_R.

Furthermore, the MoVIT algorithm was compared to straightforward extensions of FDK and SIRT. In the following, the method is demonstrated with FDK reconstructions. An equivalent strategy was also used for the SIRT reconstructions. In a first step, the conventional FDK reconstructions are registered to each other:

α_{r r^{'}} = {argmin}_{α} (\frac{1}{N} \sum_{i = 1}^{N} {(f_{r^{'}, i} - {(τ (α) f_{r})}_{i})}^{2})

where α_rr′ are the b-spline parameters describing the deformation from time frame r to time frame r′. Subsequently, the mean reconstructions f_m,1, f_m,2, . . . , f_m,R are then calculated as follows:

f_{m r^{'}} = \sum_{r \in 𝒩_{r}} z_{r r^{'}} τ (α_{r r^{'}}) f_{r}

, where

z_{r r^{'}} = exp [- {(k_{r r^{'}} / b)}^{2}] {(\sum_{l \in 𝒩_{r}} exp [- {(k_{r l} / b)}^{2}])}^{- 1}

and

k_{r r^{'}} = N^{- 1} \sum_{i = 1}^{N} {(f_{r^{'}, i} - {(τ (α_{r r^{'}}) f_{r})}_{i})}^{2}

. This technique combined with SIRT or FDK will be referred to as SIRTmean and FDKmean, respectively.

Lastly, the MoVIT method, as described in section 2, with the result of SIRTmean c_mr as initialisation will be used.

3.2. Figures of merit

The methods were evaluated using three figures of merit: mean squared error (MSE), the structural similarity index (SSIM) and feature similarity index (FSIM).

MSE: The MSE can be calculated with $N^{- 1} \sum_{n = 1}^{N} {(s_{n} - t_{n})}^{2}$ , where s is the calculated reconstruction and t is the ground truth.
SSIM: While MSE is a method which calculates absolute errors, SSIM quantifies the similarity of images as perceived by the human visual system [25]. The SSIM of a reconstruction is calculated as described in Wang et al. [25]. Since SSIM is a similarity measure, a perfect reconstruction has a SSIM of 1, the worst reconstruction a SSIM of 0.
FSIM: Similar to SSIM, FSIM is a method quantifying the perceived similarity of pictures [26]. In contrast to SSIM, it takes mainly the phase congruency and gradient magnitude into account.

3.3. Numerical simulations

Several numerical experiments were performed on a numerical phantom of a viscoelastic, open cell PU foam under compression. These models were provided by Huntsman (Everberg, Belgium) and are based on finite element simulations of different stages in the compression process. Four models were voxelised on an isotropic voxel grid of 400 × 400 × 400 from which projections of size 100 × 100 were generated, in order to avoid the inverse crime of generating the data with the same model as the one used in the reconstruction [27]. The projections were simulated in an interleaved scanning protocol with 20 projections per time frame, where each time frame has an angular range of 180 degrees [28]. Poisson distributed noise was applied on the data, assuming an incoming beam intensity of 10⁴ (photon count). During each time frame, the foam was compressed another 1.75% of the original sample height. The sample was reconstructed with a range of different reconstruction techniques as described in subsection 3.1.

The volume registration was performed with a b-spline deformation model with a control point spacing of 8 voxels. The b-spline parameters were optimized by minimizing the MSD in a multi-resolution framework. The SIRT and MoVIT algorithms were iterated until the lowest MSE was achieved. The optimal value of the parameter b, in terms of the MSE of the reconstructions, was 0.8 and was selected in this experiment. Renderings of the different reconstructions are shown in Fig. 4.

Fig. 4 Renderings of the simulated reconstruction of the compression of a foam sample with different (a) FDK, (b) SIRT, (c) FDKmean, (d) SIRTmean, and (e) MoVIT. The red circles indicate example areas where the struts that are better reconstructed in the MoVIT reconstructions in comparison with the SIRTmean reconstruction.

Download Full Size | PDF

Both the SSIM, FSIM and the MSE of the reconstructions were calculated in function of the photon count (I₀) and the number of projections per time frame (proj/time frame). These results are shown in Fig. 5 and Fig. 6, respectively.

Fig. 5 MSE, SSIM and FSIM of the reconstructions of the numerical phantoms (see Section 3.3) in function of the photon count I₀ and with 20 projections per time frame.

Download Full Size | PDF

Fig. 6 MSE, SSIM and FSIM of the reconstructions of the numerical phantoms (see Section 3.3) in function of the projections per time frame and a photon count of 10⁴.

Download Full Size | PDF

3.4. Polyurethane dataset

A dynamic X-ray CT dataset was acquired by Inside Matters with a gantry-based high-resolution scanner [29]. A viscoelastic, open cell PU sample (provided by Huntsman) of 11 mm high was loaded in a compression stage which was mounted in the scanner. Each dataset (= time frame) consisted out of 2000 equiangular projections (1316 × 1312 pixels, pixel size 0.1 mm, tube voltage 65 kV, exposure time 35 ms) acquired over an angular range of 360 degrees. Between these scans, the sample was compressed l × 0.5 mm, where l = 0, .., L − 1 is the time frame number. All reconstructions were calculated on a 1316 × 1316 × 401 isotropic voxel grid with a voxel size of 16 μm. Each time frame was reconstructed with four different methods: 1) conventional SIRT with 2000 projections/time frame, 2) conventional SIRT with 1000 projections/time frame, 3) SIRTmean with 1000 projections/time frame and, 4) MoVIT with 1000 projections/time frame. MoVIT and SIRTmean estimate the deformation and includes the projections of a single neighbouring time frame (time frame r + 1, except for the last time frame which uses time frame r − 1) to the reconstruction of a particular time frame. The SIRT reconstruction with 2000 projections has thus the best possible quality that the MoVIT reconstruction can achieve by incorporating the projections of a single neighbouring time frame. Reconstructions with only 1000 projections were performed with projections with projection numbers 1 + (r mod 2), 3 + (r mod 2), . . . , M, where r is the time frame number and M the total number of acquired projections/time frame. As such neighbouring time frames have interleaved projections. The volume registration was performed with a b-spline deformation model with a control point spacing of 8 voxels. The b-spline parameters were optimized by minimizing the MSD in a multi-resolution framework. Based on emperical evaluation the parameter b was chosen to be 0.6, which gave reasonable results. The MoVIT reconstruction was initialized with the SIRTmean reconstruction, after which 50 MoVIT iterations were performed. In Fig. 7, the horizontal cross sections of the SIRT (1000 projections/time frame), SIRTmean and MoVIT reconstructions (2 time frames) are compared with the SIRT reconstruction with 2000 projections/time frame, which serves as the ground truth, with three different metrics (MSE, SSIM and FSIM). The standard deviation of the metrics was determined by calculating the metrics on 100 horizontal cross sections. Furthermore, Fig. 7(d) shows the histogram of the region showed in Fig. 8 for the different reconstructions of the third time frame. In Fig. 8 and Fig. 9, a vertical and horizontal cross section, reconstructed with SIRT, SIRTmean, MoVIT with one additional neighbouring time frame and MoVIT with two additional neighbouring time frames (time frame r − 1 and r + 1), of the third time frame are shown, respectively. To the last two reconstructions we will respectively refer to as MoVIT (2 time frames) and MoVIT (3 time frames). In Table 1, the standard deviation of the noise at the third time frame is reported for the described methods. This standard deviation was measured selecting a region (16106 voxels) inside of a big pore.

Fig. 7 The mean MSE (a), SSIM (b) and FSIM (c) for the SIRT (yellow), SIRTmean (purple) and MoVIT with 2 time frames (green) reconstruction with 1000 proj/time frame (compared with the SIRT reconstruction with 2000 proj/time frame) for each time frame. The standard deviation of the metrics was determined by calculating the metrics on 100 horizontal cross sections. (d) shows the histogram of the region displayed in Fig. 8 for the different reconstructions of the third time frame.

Download Full Size | PDF

Table 1. Standard Deviation of the Noise

View Table | View all tables in this article

4. Discussion

4.1. Numerical simulations

Figure 4 shows renderings of the second time frame of the compressed foam sample computed with FDK (a), SIRT (b), SIRTmean (c), FDKmean (d) and MoVIT (e). FDK produces the worst results: the reconstruction is very noisy and the foam structures are very hard to differentiate from their surroundings. The FDKmean algorithm improves this reconstruction significantly, however high levels of noise are still present around the pressure plates due to cone beam artefacts. The SIRT and SIRTmean reconstructions show the foam structure nicely, however some struts of the foam are missing in the neighbourhood of the pressure plates. The MoVIT reconstruction shows the foam structure clearly and is able to reconstruct the struts close to the pressure plate better than the SIRT and SIRTmean reconstructions (see red circles in Fig. 4 for example regions). In Fig. 5, the MSE, SSIM, and FSIM are shown in function of the photon count. With respect to the MSE and SSIM, the MoVIT reconstruction performs best. This is also the case for the FSIM metric at high (> 10⁴) photon counts or a high (> 30) number of projections per time frame. The SIRTmean method provides slightly worse reconstruction in terms of MSE and SSIM, but is an improvement over the conventional SIRT reconstructions. For low photon counts and low number of projections per time point, the FSIM of the SIRTmean are slightly better than the MoVIT reconstruction. Only at very low noise levels the MoVIT method results in a slightly worse reconstruction quality compared to the SIRT methods, which is caused by inaccurate deformation estimations on the very low SNR initial SIRT reconstructions. From Fig. 6, similar conclusions can be drawn. Here the MoVIT method provides the best results with respect to the MSE and SSIM.

4.2. Polyurethane dataset

The results of the polyurethane foam dataset are shown in Fig. 8 and Fig. 9. The metrics of the different reconstructions in comparison with the SIRT reconstruction with 2000 proj/time frame are shown in Fig. 7. The SIRT reconstruction with 2000 projections/time frame shows the foam structures very clearly. The SNR is sufficiently high to observe small foam structures. As can be expected, the noise of the SIRT reconstruction with only 1000 projections/time frame has a higher standard deviation than the SIRT reconstruction with 2000 projections/time frame and image details are obscured by the noise. The noise in the SIRTmean reconstruction (1000 projections/time frame) has a larger standard deviation as the SIRT reconstruction with 2000 projections/time frame, however it is considerably lower than that of the SIRT reconstruction with 1000 projections/time frame. The MoVIT method is able to lower the standard deviation of the noise even further. This is also reflected by the histograms in Fig. 7(d), which show a narrower background peak for the methods with a smaller standard deviation. By taking an extra time frame into account the results of the MoVIT algorithm improve even further. The standard deviation of the noise is even lower of that of the SIRT reconstruction with 2000 projections/time frame. Additionally the MSE, SSIM and FSIM metrics (see Fig. 7) reveal that the MoVIT reconstruction is more similar to the SIRT reconstruction with 2000 projections/time frame than the SIRTmean reconstruction. These results indicate that the improvement of MoVIT algorithm is not simply an averaging effect which improves the reconstruction. Indeed, incorporating more projections into the MoVIT reconstruction process not only improves the SNR but also improves the spatial resolution.

Fig. 8 Zoomed horizontal cross section through the third time frame of the polyurethane dataset.

Download Full Size | PDF

Fig. 9 Vertical cross section through the third time frame of the polyurethane dataset.

Download Full Size | PDF

5. Conclusion

In this paper, we have presented the MoVIT frame work which aims to reconstruct dynamic CT datasets with high temporal resolution and spatial resolution. The MoVIT frame work estimates the deformation between different time frames and enables including the projections of these time frames in the reconstruction of a certain time frame without introducing deformation artefacts. The method was validated on numerical simulations and a real dataset of polyurethane foam under compression. Both experiments showed an increase of reconstruction quality with respect to conventional reconstruction algorithms. The performed experiments show that the MoVIT algorithm is able to successfully exploit the data redundancy present in 4D-CT datasets. It allows lowering the acquisition time of a single time frame without compromising the SNR of the reconstructed images.

Funding

Fund for Scientific Research-Flanders (FWO) (S004217N, G0F9117N); Agency for Innovation by Science and Technology in Flanders (IWT)(SBO 120033); iMinds-ICON (130647).

Acknowledgments

The authors would like to thank Maarten Moesen and Mark Brennan (Huntsman) for providing the numeric polyurethane models and Jelle Vlassenbroeck (Inside Matters) and the Centre for X-ray Tomography of the Ghent University (UGCT) for the acquisition of the polyurethane datasets.

References and links

1. J. Elliott, A. Windle, and J. Hobdell, “In-situ deformation of an open-cell flexible polyurethane foam characterised by 3D computed microtomography,” J. Mater. Sci. 37, 1547–1555 (2002). [CrossRef]

2. A. T. Jang, J. D. Lin, Y. Seo, S. Etchin, A. Merkle, K. Fahey, and S. P. Ho, “In situ compressive loading and correlative noninvasive imaging of the bone-periodontal ligament-tooth fibrous joint,” J. Vis. Exp. 85, 51147 (2014).

3. F. Koksel, S. Aritan, A. Strybulevych, J. H. Page, and M. G. Scanlon, “The bubble size distribution and its evolution in non-yeasted wheat flour doughs investigated by synchrotron X-ray microtomography,” Food Res. Int. 80, 12–18 (2016). [CrossRef]

4. G. Van Eyndhoven, K. J. Batenburg, D. Kazantsev, V. Van Nieuwenhove, P. D. Lee, K. J. Dobson, and J. Sijbers, “An iterative CT reconstruction algorithm for fast fluid flow imaging,” IEEE Trans. Image Process. 24, 4446–4458 (2015). [CrossRef] [PubMed]

5. M. Negahdar and A. A. Amini, “Estimation of affine motion from projection data using a mass conservation model,” in Proceedings of Annual International Conference of the IEEE Engineering in Medicine and Biology Society (IEEE, 2011), pp. 8041–8044.

6. W. Wein and A. Ladikos, “Towards general motion recovery in cone-beam computed tomography,” in Proceedings of Fully Three-Dimensional Image Reconstruction Radiology and Nuclear Medicine (2013), pp. 54–57.

7. R. Frysch and G. Rose, “Rigid motion compensation in interventional C-arm CT using consistency measure on projection data,” in Proceedings International Conference on Medical Image Computing and Computer-Assisted Intervention (Springer International Publishing, 2015), pp. 298–306.

8. V. Van Nieuwenhove, J. De Beenhouwer, T. De Schryver, L. Van Hoorebeke, and J. Sijbers, “Data-driven affine deformation estimation and correction in cone beam computed tomography,” IEEE Trans. Image Process. 26, 1441–1451 (2017). [CrossRef] [PubMed]

9. G. R. Myers, A. M. Kingston, T. K. Varslot, M. L. Turner, and A. P. Sheppard, “Dynamic tomography with a priori information,” Appl. Opt. 50, 2685–3690 (2011). [CrossRef]

10. G. Van Eyndhoven, K. J. Batenburg, and J. Sijbers, “Region-based iterative reconstruction of structurally changing objects in CT,” IEEE Trans. Image Process. 23, 909–919 (2014). [CrossRef] [PubMed]

11. D. Kazantsev, W. M. Thompson, W. R. B. Lionheart, G. Van Eyndhoven, A. P. Kaestner, K. J. Dobson, P. J. Withers, and P. D. Lee, “4D-CT reconstruction with unified spatial-temporal patch-based regularization,” Inverse Probl. Imag. 9, 447–467 (2015). [CrossRef]

12. K. A. Mohan, S. V. Venkatakrishnan, J. W. Gibbs, E. B. Gulsoy, X. Xiao, M. De Graef, P. W. Voorhees, and C. A. Bouman, “TIMBIR: a method for time-space reconstruction from interlaced view,” IEEE Trans. Comp. Imaging 1, 96–111 (2015). [CrossRef]

13. M. Brehm, P. Paysan, M. Oelhafen, P. Kunz, and M. Kachelrieß, “Self-adapting cyclic registration for motion-compensated cone-beam CT in image-guided radiation therapy,” Med. Phys. 39, 7603–7618 (2012). [CrossRef] [PubMed]

14. A. A. Isola, M. Grass, and W. J. Niessen, “Fully automatic nonrigid registration-based local motion estimation for motion-corrected iterative cardiac CT reconstruction,” Med. Phys. 37, 1093–1109 (2010). [CrossRef] [PubMed]

15. H. Yan, X. Zhen, M. Folkerts, Y. Li, T. Pan, L. Cervino, S. B. Jiang, and X. Jia, “A hybrid reconstruction algorithm for fast and accurate 4D cone-beam CT imaging,” Med. Phys. 41, 071903 (2014). [CrossRef] [PubMed]

16. C. P. V. Christoffersen, D. Hansen, P. Poulsen, and T. S. Sorensen, “Registration-based reconstruction of four-dimensional cone beam computed tomography,” IEEE Trans. Med. Imaging 32, 2064–2077 (2013). [CrossRef] [PubMed]

17. J. Gregor and T. Benson, “Computational analysis and improvement of SIRT,” IEEE Trans. Med. Imaging 27, 918–924 (2008). [CrossRef] [PubMed]

18. D. Rueckert, L. I. Sonoda, C. Hayes, D. L. G. Hill, M. O. Leach, and D. J. Hawkes, “Nonrigid registration using free-form deformations: application to breast MR images,” IEEE Trans. Med. Imaging 18, 712–721 (1999). [CrossRef] [PubMed]

19. M. Chen, W. Lu, Q. Chen, K. J. Ruchala, and G. H. Olivera, “A simple fixed-point approach to invert a deformation field,” Med. Phys. 35, 81–88 (2008). [CrossRef] [PubMed]

20. B. K. Bay, T. S. Smith, D. P. Fyhrie, and M. Saad, “Digital volume correlation: Three-dimensional strain mapping using X-ray tomography,” Exp. Mech. 39, 217–226 (1999). [CrossRef]

21. W. J. Palenstijn, K. J. Batenburg, and J. Sijbers, “Performance improvements for iterative electron tomography reconstruction using graphics essing units (GPUs),” J. Struct. Biol. 176, 250–253 (2011). [CrossRef] [PubMed]

22. W. Van Aarle, W. J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt, A. Dabravolski, J. De Beenhouwer, and J. Sijbers, “Fast and flexible X-ray tomography using the ASTRA toolbox,” Opt. Express 24, 25129–25147 (2016). [CrossRef] [PubMed]

23. W. J. Palenstijn, J. Bédorf, and K. J. Batenburg, “A distributed SIRT implementation for the ASTRA toolbox,” in Proceedings Fully Three-Dimensional Image Reconstruction Radiology and Nuclear Medicine (2015), pp. 166–169.

24. L. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical cone-beam algorithm,” J. Opt. Soc. Am. A 1, 612–619 (1984). [CrossRef]

25. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Trans. Image Process. 13, 600–612 (2004). [CrossRef] [PubMed]

26. L. Zhang, L. Zhang, X. Mou, and D. Zhang, “FSIM : A feature similarity index for image quality assessment,” IEEE Trans. Image Process. 20, 2378–2386 (2011). [CrossRef] [PubMed]

27. J. Kaipio and E. Somersalo, “Statistical inverse problems: Discretization, model reduction and inverse crime,” J. Comput. Appl. Math. 198, 493–504 (2007). [CrossRef]

28. A. Kaestner, B. Münch, P. Trtik, and L. Butler, “Spatiotemporal computed tomography of dynamic processes,” Opt. Eng. 50, 123201 (2011). [CrossRef]

29. M. Dierick, D. Van Loo, B. Masschaele, J. Van den Bulcke, J. Van Acker, V. Cnudde, and L. Van Hoorebeke, “Recent micro-CT scanner developments at UGCT,” Nucl. Instr. Meth. Phys. Res. B 324, 35–40 (2014). [CrossRef]

MoVIT: a tomographic reconstruction framework for 4D-CT

Abstract

1. Introduction

2. Method

2.1. Stationary algebraic tomography

2.2. Dynamic tomographic model

2.3. Motion vector-based iterative technique

2.4. Deformation estimation

3. Experiments and results

3.1. Reconstruction methods

3.2. Figures of merit

3.3. Numerical simulations

3.4. Polyurethane dataset

4. Discussion

4.1. Numerical simulations

4.2. Polyurethane dataset

5. Conclusion

Funding

Acknowledgments

References and links

Cited By

Figures (9)

Tables (2)

Equations (10)

Optics Express

method	proj/time frame	noise standard deviation
SIRT	2000	1, 22 × 10⁻⁵
SIRT	1000	2, 41 × 10⁻⁵
SIRTmean (2 time frames)	1000	1, 63 × 10⁻⁵
MoVIT (2 time frames)	1000	1, 47 × 10⁻⁵
MoVIT (3 time frames)	1000	1, 19 × 10⁻⁵