1.

INTRODUCTION

Computed tomography (CT), as a noninvasive method, is widely used in many applications to reveal the inner structure of a target. In biomedical applications, however, the target is usually non-stationary so that the projection measurements are time-dependent. That is, the target changes during the scanning, a case commonly referred as dynamic tomography. Such changes can be periodic (e.g., the beating of a heart) or aperiodic (e.g., the flow of contrast agent in a blood vessel), in each case, the reconstruction problem becomes severe ill-posed. If traditional static CT reconstruction methods (such as FBP, ART, IART, etc.) are directly applied, the reconstructions will be degraded seriously (e.g., being blurred or getting other artifacts), even if with small motions. A number of approaches have been proposed for dynamic CT reconstructions. For very slow or periodic movements, e.g. heart and lung imaging, gating techniques can be used.^{1, 9} That is, the scanning is executed at specific time during the periodic motion, or the acquired projection data are selected according to the periodicity. For general motions, motion compensation techniques can be used to improve the reconstruction quality.^{2–5, 7} By motion compensation, prior information about the underlying dynamic process can be utilized by the CT reconstruction algorithms. However, this technique usually assumes simple or known motion, which is not the case for a wide variety of complex motions.

Dynamic CT reconstruction aims to reconstruct a series of image slices, while each image slice undergoes very limited resources, i.e. undersampled scanning. To reduce artifacts and improve reconstruction quality, many methods have been proposed in the literature by incorporating various priors. One of the main approach is to cast the reconstruction problem under the Bayesian framework. Then, the priors could be naturally encoded by the prior distribution, and by employing the Kalman filter, priors concealed between image slices could also be easily modelled and brought into the reconstruction process. A serious drawback of this approach is that the precision matrices, which are used to encode the needed priors, are usually non-sparse huge ones. Besides, the inversion of these large-scale matrices are needed during the reconstruction process, which is not applicable for real applications. To address this issue, the dimension reduction Kalman filter (DR-KF) was introduced in⁶ to reduce the computational complexity. This is actually an effective method to serve the purpose. However, images reconstructed by this method are often smeared by “ringing” artifacts. In our experiments, even a small portion of dimension reduction leads to clear sharp ringing.

To combat the ringing artifacts and reduce the computational complexities, in this paper, we will introduce a splitting technique to approximate the priors encoded by the precision matrices with a spatial prior plus a temporal prior. This approximation helps to greatly reduce the computational complexity and totally avoid the ringing artifacts. Besides, to further improve the reconstruction quality, we propose a block-wise Kalman filter for the backward smoothing pass, which calls for much less computations and achieves very close performance compared to the original Kalman filter.

The remainder of this paper is organized as follows. In Section 2, we will briefly introduce the DRKF method for dynamic X-Ray tomography and analyze its ”ringing” artifacts. Then, the proposed method shall be described in details. Numerical experiments are performed in Section 3 to validate the proposed method. Finally, we conclude our paper in Section 4.

2.

MATERIALS AND METHODS

2.1

Dynamic X-Ray Tomography

In dynamic X-ray tomography, the attenuation variable is a function of the “time” k ∈ N. In a discrete setting, is a two-dimensional image, which is related to the observation by the following model

In tomography, is the measurement vector called sinogram, is the vectorized image, A_k is the measurement matrix and represents noise distribution. Thus, in dynamic X-ray tomography, the aim is to reconstruct a set of images changing over time.

2.2

Dimension Reduction Kalman filtering

In Bayesian framework, prior probability distributions are introduced to bring into prior knowledge. Suppose , then the posterior density can be written as⁸

In linear Kalman filtering, a linear operator is introduced to model the state transfer process. The KF is generally divided into two steps: prediction and update, which could be described as

where the matrix M_k moves the previous state to . Take the transferred state as a prediction, the complete process can be written as

Where

Direct application of the the above procedure is very time consuming. Suppose the image slices are size of N×N. The precision matrix , which are size of N²×N² dense matrices, are involved in the computations, which is not affordable even for small images like N =128. To reduce the computational complexity, the dimension reduction Kalman filter is introduced for dynamic tomography in.⁸ It starts by parameterizing the state , where P_k,r is constructed by the first r largest singular values of and the corresponding singular vectors. By setting a small r, the computational complexity shall be sharply reduced since the matrices involved in the computations are now size of r² × r². However, the SVD decompositions, which operate on the original matrices are required for each state updating, which is still very time consuming. This is solved by assuming that all variations live in the same fixed subspace, i.e. fixing P_k,r = P_0,r ≜ P_r, which can be built by performing SVD decomposition once on the initial covariance matrix . Let ∑ = USU^T, then P_r = , and (2) and (3) can be rewritten as:

where

It has been mentioned in⁶ that the reconstructed images from the dimension reduction KF might suffers from “ringing” artifacts. This is demonstrated in the following experiment.

Let the image slices are size of 128×128, and assume the standard prior Gaussian covariance matrix,⁶ which is then size of 16384 × 16384 with the (i, j)th element defined as

where σ² and l are configuration parameters, and d (x_i, x_j) denotes the Euclidean distance between pixels x_i and x_j.

Both observation and model error covariance matrices are set to σ²I and M_k = I, σ = 0.1, l = 1.5. The phantom3d function in MATLAB is employed to construct the dynamic images with 33 slices. The open-source library ASTRA is utilized to perform the forward and backward projections. There are only 5 projection views at each moment. The reconstruction results from the dimension reduction Kalman filter (DR-KF) with r = 1000 are shown in Fig.1. As shown in the last two columns, heavy ringing artifacts presents for the time-frame t =15 and 20. At the early stage, the image slices suffer from severe undersampling artifacts and the ringing artifacts are buried inside. When the image slices getting better reconstructed, the ringing artifacts also show up and then persist in subsequent reconstructions.

Figure 1.

Illustration of ringing artifacts from DR-KF. Slices of phantom3d reconstruction at time step t = 1, t = 5, t = 15, t = 20.Dimension reduction (r = 3000) methods with KF (5 angles).

To further verify the ringing artifacts, another experiment is performed to demonstrate the influence of the ratio of dimension reduction. As shown in Figure 2, by increasing the dimension parameter r, the ringing artifacts are suppressed. However, even for r = 12000, which amounts to a reduction ratio 12000/16384 ≈ 73%, the ringing artifacts are still observable. For r = 5000, which amounts to a ratio 3/10, heavy ringings present. Clearly, such a reduction ratio is far from satisfactory, since the computational complexity is still very high.

Figure 2.

The influence of the ratio of dimension reduction on ringing artifacts. Phantom3d reconstruction at time step t = 15, σ = 0.1, l = 1.5. Dimension reduction methods with r = 1000, 5000, 8000, 12000, KF(5 angles).

One may argue that the ringing artifacts could be introduced by improper choice of the configuration parameters for the Gaussian prior. So, another experiment is performed with various parameter settings. As shown in Figure 3, the parameters σ and l do have some influence on the ringing artifacts. However, they just have certain influence on the ringing patterns rather than on the magnitude of the ringings.

Figure 3.

The influence of configuration parameters of Gaussian prior on ringing artifacts. Slices of phantom3d reconstruction at time step t = 1, t = 5t = 15, t = 20. Dimension reduction methods with r = 1000 KF(5 angles).

2.3

Our Method

2.3.1

Decomposition of the prior

The computational complexity can be attributed to the regularization term

which is used to encode the prior information. The complexity of the covariance matrix results in a complex regularizer. This motivates us to approximate the above term with simpler, easier handled terms. The function Φ could always be decomposed as

The first two terms can be though of being encoding the priors regarding to and , respectively, while the last term expresses the couplings between adjacent image slices. Since the prior regarding to is not relevant for the reconstruction of , the second term could be removed. So, we propose the following posterior distribution

The regularization term represents the prior knowledge regarding to the current ideal image . This prior is usually available, e.g. CT images are usually assumed to be piecewise constant and thus possess sparse gradients. The regularization term models the continuity along the temporal axis. A simple strategy is to just take as a prior reference image.

So, we are proposing to decompose the regularization term involving the covariance matrix into a spatial regularization term plus a temporal regularization term. By choosing these two regularizers properly, the computations could be greatly reduced and ringing artifacts could also be completely avoided.

2.3.2

The online forward pass

By choosing total variation regularization for the spatial regularizer, we propose the following model for dynamic tomography reconstruction

where the parameter λ weights the importance of the prior image while β weights the spatial prior. To solve the above minimization problem, the primal-dual based Chambolle-Pock (CP) algorithm is employed.

2.3.3

Block Kalman Filter

In the Kalman filtering framework, a backward smoothing procedure is usually applied for offline reconstructions. When all the projection data have been acquired and all image slices have been reconstructed, then the Kalman filter could be applied once more reversely starting from the last image slice. The backward smoothing procedure could significantly further improve the reconstruction quality.

However, the backward filtering suffers from the same drawbacks as the forward filtering. To take advantage of the backward smoothing, we propose two approximations to reduce the computational complexity. The first one is to replace the system matrix A_k by identity matrices. This is reasonable since all the image slices have been reconstructed. The second one is to decompose the filtering into blocks, i.e. decomposing the image slices into overlapping blocks and applying the Kalman filter on each block. The final results are synthesized from the filtered blocks.

With much less computations, the block Kalman filter could attain almost the same smoothing results as the non-block version does. For example, for an image size of 256 × 256, we can divide it into 81 overlapping blocks such that each block is size of 32×32. In the block Kalman filtering process, the covariance matrix size of 65536×65536 is replaced by 81 image blocks size of 1024×1024. Clearly, storage and computational complexities shall be dropped considerably.

3.

NUMERICAL EXPERIMENTS

Similarly, in the phantom3d simulated data, which is size of 256×256×256, we take 50 consecutive slices in the vertical direction. The projection data are acquired with 10 equally distributed scanning angles for each slice.

We set R_k = I,Q_k = 0.1 * I. The experiments are carried out on a Linux server with two Intel(R) Xeon(R) Gold 6132 CPUs @ 2.60GHz and 256G memory. The server is also equipped with 8 Nvidia GeForce RTX 2080 Ti GPU cards, and ASTRA (http://www.astra-toolbox.com/) uses one of them for applying the forward and backward projectors.

The reconstruction results are illustrated in Fig. 4. As shown in the first row, the ringing artifacts have been completely avoided in the image slices reconstructed by the proposed forward pass, i.e. the online reconstruction procedure. The backward smoothing could help to reduce noise, as shown in the second row. However, structuring information could also be wrongly propagated back during the backward smoothing pass. In fact, we also implemented the non-block version of the Kalman filtering for the backward smoothing pass, and the results of block version is a little better than those of the non-block. This can be told from the quantitative indices listed in Table 1. It should be mentioned that, to process one image, the block-KF takes 7s for one image, while the non-block KF takes about 4000s.

Figure 4.

Different reconstructions of the shepp-logan phantom with the forward pass(top), the block-KF smoothing(middle), and the reference imagse at time step t = 1,10,30.

Table 1.

PSNR of different methods for different frames

4.

CONCLUSION

To reduce the computational complexity of the Kalman filtering approach for dynamic X-ray tomography, we propose to approximate the regularization term involving the covariance matrix with two simpler regularization terms, one accounts for spatial correlation while the other one accounts for temporal correlations. The resulting method calls for much less storage and computational complexity while achieves competitive reconstructions. Besides, a block version of the Kalman filter is proposed for the backward smoothing procedure, which achieves similar smoothing behavior compared to the non-block one while again requires much less computational resources.

In our numerical experiments, the matrix M_k is always set as the identity matrix. In fact, we can also use the optical flow method to compute a better M_k such that the reconstruction quality could be further improved. This shall be investigated soon.

In order to further improve the computation speed, we also plan to use GPU and deep learning methods to approximate the forward reconstruction pass.