Temporal sparsity exploiting nonlocal regularization for 4D computed tomography reconstruction

2.1Nonlocal regularization on weighted graphs for time-lapse tomography

By taking a nonlocal regularization approach based on weighted graphs [21] we can rewrite the second minimization problem in (3) as:

∇ _ωX (v) ∥ is the weighted local variation of X over the graph. The ω-weighted p-Laplace term in (4) is convenient to operate with since various choices of p lead to well-known nonlocal filters used in image processing [22, 23]. For instance, p = 2 leads to the weighted linear heat diffusion on graphs and p = 1 to the mean curvature or Total Variation (TV) seminorm on graphs [21].

Two vertices u and v are said to be adjacent if the edge (u, v) ∈ E connects them, the weight between u and v is denoted by ω (u, v): ∀ (u, v) ∈ E, ω (u, v) = ω (v, u). Also ω (u, v) >0 if u ≠ v and 0 otherwise. Notation u ∼ v will be used for two adjacent vertices. The vertex u ∈ B (v), where B (v) is a cubic box (the searching space) size of s_x × s_y × s_k. Potentially, the box size can be extended to the 4D case (s_x × s_y × s_z × s_k), however due to significant computational costs involved it is currently infeasible for large XMT data. Therefore, our regularization is currently in 3D (x, y+time) and it has similarities with the video denoising problem [26]. The major benefit of this realization is that it is relatively fast and can be programmed in a “slice-by-slice” manner which does not require lots of computer memory for reconstruction and is easily parallelizable.

For regularization on graphs we define nonlocal derivative operators. The local variation of the weighted gradient operator

∇ _ωX (v) ∥ is defined by

The weighted p-Laplace operator at a vertex v is defined as

One can see that for the case p = 2: γ (u, v) = ω (u, v) and for p = 1 it is:

Both cases are convex and therefore lead to the convex minimization problem (4), which guarantees a unique solution X˜ under a proper choice of parameter β. In this paper we use p = 1 which provides a better regularization performance than the linear model p = 2. Note that the nonlocal TV seminorm (5) is smoothed by a small constant to ensure differentiability: ∥∇ωX(v)∥ε=∑u∼v(∂uX(v))2+ɛ2. The non-convex case p < 1 is not considered in this paper [21].

Therefore for p ≥ 1 the system (4) has a unique solution written as:

The linearized Gauss-Jacobi iterative method can be used to solve the system (10) which also can be considered as the discrete Euler-Lagrange equation of minimization problem (4) [21–23]. Let t be an iteration step, then the fixed point (FP) iteration is used to solve (10):

In our case (p = 1), γ (u, v) as taken as in (8). The iteration procedure (11) can be stopped using criterion X^t+1 -  X^t ∥ < ς, where ς is a small constant. In our experiments we use only one iteration of (11) mainly due to regularization performance restrictions.

The regularization on graphs for time-lapse reconstruction involves a search for similar vertices in a box B (v). Similar to nonlocal video denoising [26], a quadratic patch F size of r_x × r_y is used to provide more robust estimation of similarity between vertices u and v (calculation of Euclidian distance on patches [22, 23]).

2.2Temporal sparsity exploiting nonlocal regularization for 4D tomography

The computational complexity of the proposed algorithm can go up to O(N²) and therefore becomes infeasible for large XMT datasets. Here we propose a simple approach to significantly accelerate the regularization process on graphs by considering sparsity in the temporal domain. We developed an accelerated regularization on graphs (ARG) technique based on two observations: a) the temporal data are likely to be sparse and redundant (some features remain stationary or change only a little over time); b) two vertices are likely to be similar (structurally) if the difference in intensity between them is lower than the maximum level of noise present in the signal. We use these two observations to modify and accelerate the classical regularization algorithm on graphs (RG) [21].

First we consider our (a) observation of the temporal sparsity. We say that the s_x,y size of the searching box B (v) must be larger if a feature at the vertex v is changing (e.g. moving) quickly over time and s_x,y of B (v) must be smaller if the feature is almost stationary (e.g. a uniform region or an edge). By reducing B (v) for temporally non-active (stationary) features we drastically reduce our search space. In order to decide if the feature is stationary or not and if the vertex u is similar to v we use the observation (b).

While minimizing the cost function (2), we intentionally do not impose any positivity constraints on our solution and therefore Xⁿ will have negative values on every n-iteration. The maximum noise level present in our reconstructed image can be simply found as δⁿ = | min( Xⁿ) |. We also say that the trustworthy intensity differences between vertices u and v are bounded above by Lδⁿ, where L ∈ (0, 1] is an empirical constant (we use L = 0.4 for all our experiments). To establish the optimal searching box size B_o (v) for v (i, j, k) ,  k = 1, 2, …, K we calculate the following measure of temporal sparsity:

Then, if we select the lower size of the box B_o (v) as sxl=syl and the upper size as sxu=syu, S (v) can be linked to B_o (v) (see Fig. 1 (right)). We divide the whole interval min(S_k) - max(S_k) into an equal number of D steps (we use D = 10 in our experiments) and linearly assign values for B_o (k) in a way that min(S_k) corresponds to s^u and max(S_k) to s^l respectively. The simple reasoning here is that the larger S values should be related to the smaller sizes of the searching window and vice versa.

Finally, the following condition is used to calculate nonlocal weights in (11):

One can see that we reduce the data space substantially by first creating a spatially variant search box B_o and next we remove vertices which do not fit our assumptions of similarity using the rule (13). As we will demonstrate later, this approach provides an acceleration by an order of magnitude without impairment of reconstruction quality. Moreover, the averaging across dissimilar features in ARG is much reduced and leads to improved resolution.

3.1Tomographic reconstruction (synthetic data)

To test ARG and RG approaches we created a dynamically evolving phantom from high quality tomographic measurements (the phantom is available here [27]). The high resolution data of a mouse tibia sample was acquired at the Diamond-Manchester branchline (I13) of the Diamond Light Source (DLS). In order to simulate the dynamic behaviour of the stationary sample we take 10 vertical slices of 3D reconstructed volume and shift them horizontally using an equidistant step of 15 pixels (see Fig. 2).

Three time frames (k = 1, 5, 10) of the mouse tibia phantom are shown in Fig. 2. The total size of the phantom is 400 × 400 × 10 (K) pixels. To avoid the “inverse crime” of reconstructing on the same grid where projection data was simulated, we used a higher resolution of the phantom on a 800 × 800 isotropic pixel grid to generate projections with a strip kernel [24]. Gaussian noise was applied to the projection data (5% of a true signal maximum). Reconstructions were calculated on a smaller 400 × 400 isotropic pixel grid and with a linear projection model [24]. We used 180 and 90 projection angles in [0, π) angular interval (assuming parallel beam geometry) to reconstruct the phantom. The Root Mean Square Error (RMSE) is calculated as

In Fig. 4 we show reconstructions (for time frame k = 5 of the phantom in Fig. 2 (middle)) with various methods for two cases of limited noisy projections (180 and 90 projections). The main aim here is to show that the accelerated ARG method can be competitive (quantitatively) with the classical RG method, while much faster computationally. To reconstruct the data we used various sizes (s_x = s_y) of the searching window B_sx,sy,sk, while s_k = 9 remains constant (we used nine temporal frames for every time-frame in regularizion). Because our dynamic phantom (see Fig. 2) has significant vertical shifts in time, one should use larger search window sizes (s_x = s_y) to track the similar edge. The RMSE was calculated for all reconstructed phantoms and given in Tables 2 and 3.

One can see that the reconstructions with the CGLS method are noisy (see Fig. 4) and have the highest errors. The regularization on graphs (RG) is able to reduce noise significantly and sharpen important features. However, the small search window reconstructions (method RG-B_{9, 9, 9}) have a distinctive blocky appearance (insufficient amount of temporal information collected). With a larger search window (RG-B_{43, 43, 9}), reconstructions are much smoother and edges are well preserved. Unfortunately, the computational time is very large for RG-B_{43, 43, 9} (see Tables 2, 3) and therefore the method is impractical. The computation time is drastically reduced with the proposed method (ARG-B_o,9) while the error is the lowest. Notably, the quality of reconstruction is very similar to the original method (RG-B_{43, 43, 9}), however some small and thin features are better preserved than with RG-B_{43, 43, 9}. The effect of under and over estimation for various methods can be clearly seen in a 1D plot shown in Fig. 5 (by arrows we show the most erroneous regions in the profile recovered by methods RG-B_{9, 9, 9} and RG-B_{43, 43, 9}). Notably, the RG-B_{43, 43, 9} method smooths small features while the proposed method recovers them much better (see the parabola shaped region in Fig. 5).

From those numerical experiments we conclude that the proposed method (ARG-B_o) is more feasible for larger datasets due to performance and quality-wise it is similar or better than the classical nonlocal regularization approaches.

3.2Reconstruction of a real XMT dynamic dataset

The main purpose of developing such a spatio-temporal regularization technique (see Section section 2) is to apply it to time-lapse XMT data. Due to high resolution detectors and rapid acquisition the data sizes are very large [10, 11]. Therefore fast numerical implementation of the projection operator A and regularizer R ( X) is very important for iterative minimization of (2). In order to perform fast projection-backprojection operations we use GPU-accelerated open-source package: the ASTRA toolbox [28]. The proposed regularizer is implemented in C [27], however it can be also realized on a GPU which can provide substantial acceleration.

The following time-lapse experiment was performed on the I13 beamline of DLS, where multiple fast 3D scans of time-varying sample (melting-freezing ice-cream) were obtained. The ice-cream microstructure consists of three phases, air cells (gas), ice crystals (solid) and unfrozen matrix (liquid) (see Fig. 6, 7). The aim of the experiment was to understand how phase content varies with temperature variation [5]. During the thermal cycling the inner structure of a sample is changing (crystals melt or coarsen). Since dynamic changes are easily traceable (controlled sample cooling and heating procedure), this dataset is a perfect case study for our reconstruction technique.

Additionally, since our reconstruction method searches for some structural correlations in time it significantly enhances the edges between different phases. Due to poor contrast between ice-crystals and air bubbles one has to recover the edges of the unfrozen matrix as accurately as possible.

In Fig. 6 one can see reconstructions with three different methods: FBP (Filtered Back-Projection [24]), CGLS and the proposed ARG method. The time-lapse data size is 2k × 2k × 1k × 6 (K) voxels, and therefore 6 large 3D volumes are reconstructed. One axial slice from the first volume is presented in Fig. 6.

To reconstruct the data iteratively we perform 40 iterations of CGLS method with one FP iteration (11) for the ARG penalty. All parameters for the ARG method except β and h are selected the same as in Table 1. The computation time of one ARG iteration (11) to regularize 6 slices (2k × 2k × 6 (K)) is 170 seconds (all available time frames are used). The classical RG implementation is not used here due to long computation times (more than 1500 seconds for one iteration).

From the reconstructions it can be seen that the direct method (FBP) and CGLS both deliver unquantifiable reconstructions due to high noise levels and ring artifacts present. It can be seen in Fig. 6 that the ARG method is able to suppress noise significantly while also enhance edges between different phases. This enhancement is critical for the subsequent segmentation (see Fig. 7), where different phases can be extracted and quantified.