Deep Inertia $L_{p}$ Half-quadratic Splitting Unrolling Network for Sparse View CT Reconstruction

In order to take advantage of the sparsity of the CT image, we use the $L_{p}$-norm with the wavelet transform to form the regularization term in (2), i.e.

where $0<p<1$, $W=(W_{1},W_{2},...,W_{m})$ is a $m$-channel operator with $W_{i}^{T}W_{i}=I(1\leq i\leq m)$. The operator $W$ is chosen as the highpass components of the piecewise linear tight wavelet frame transform [28]. Solving (3) directly is challenging due to the discontinuity in the $L_{p}$ norm. To solve (3), we introduce the auxiliary variable $z$, and use the half-quadratic splitting method to obtain the following problem:

where $\lambda>0$ and $\gamma=(\gamma_{1},\gamma_{2},...,\gamma_{m})$ with $\gamma_{i}>0$ $(i=1,2,...,m)$ are weight parameters. The above problem can be split into $z$-subproblem and $u$-subproblem solved in alternating iterations. To speed up the convergence, inspired by the inertia algorithm [29, 30], we introduce inertia steps in solving the two subproblems, i.e.,

where (5) and (7) are subproblem solving steps, (6) and (8) are inertia steps, and $\alpha\in[0,1)$ and $\beta\in[0,\frac{\sqrt{5}-1}{2})$ are inertia weight parameters.

The $z^{k+1}$-update step (5) can be solved by the proximal operator ${\rm prox}_{p,\eta}$, which is defined as

where $0<p<1$, it can be computed as [31]

where $\rho=(2(1-p)/\eta)^{\frac{1}{2-p}}$, $\tau=\rho+p\rho^{p-1}/\eta$, $g$ is the solution of $h(g)=pg^{p-1}+\eta g-\eta|t|=0$ over the region $(\rho,|t|)$. Since $h(g)$ is convex, when $|t|>\tau$, $g$ can be efficiently solved using Newton’s method.

For the $u^{k+1}$-update step (7), according to the Euler-Lagrange equations, we have:

To avoid the inverse, this equation can be solved using the conjugate gradient method. We call our algorithm as Inertial HQS_p-CG (IHQS_p-CG) for short and summarized in Algorithm 1. Next, we establish the global convergence of the IHQS_p-CG algorithm.

In this section, we provide a convergence analysis of Algorithm 1.

Details of the proof can be found in the supplementary material.

The deep learning algorithm used training and validation data from the “2016 NIH-AAPM-Mayo Clinic Low Dose CT Grand Challenge” dataset [35]. The training set had 1596 slices from five patients, and the validation set had 287 slices from one patient. Experiments were uniformly performed using a fan beam geometry with 800 detector elements. The test data included 38 slices from new patients in the AAPM dataset and 30 slices from the Covid-19 dataset [36] to assess the algorithm’s generalization and robustness.

Our proposed algorithms were compared with eight others, including five traditional (FBP [3], BM3D [5], HQS-CG [13], combining compressed sensing and sparse convolutional coding for PWLS-CSCGR [10], HQS_p-CG without inertial $L_{p}$ norms) and three deep learning algorithms (PDNet [17], FBPNet [15], MetaInvNet [13]). For our lightweight IHQS_p-Net* and IHQS_p-Net, we set K = 6 and limited the CG algorithm iterations to 7. Data were degraded with different angles (60, 90, 120, 180) and two noise levels (Gaussian noise $\sigma=0.3$, mixed Gaussian and Poisson noise with intensity $I_{0}=5\times 10^{5}$).

Our models were implemented using the PyTorch framework. The experiments were conducted on PyTorch 1.3.1 backend with an NVIDIA Tesla V100S GPU, 8-core CPU and 32GB RAM. For IHQS_p-CG, the hyperparameters were set to $p=0.7,\alpha=0.5$, $\beta=0.6$, $\lambda\in(0.0004,0.0008)$, $\gamma\in(0.1,0.4)$. For IHQS_p-Net, the hyperparameters were set to $p=0.7,\alpha=0.3$, $\beta=0.1$, $\lambda\in(0.001,0.004)$, $\gamma\in(0.01,0.06)$.

Evaluation metrics used are peak signal-to-noise ratio (PSNR) [37] and structural similarity (SSIM) [34]. For a comprehensive evaluation of the algorithm performance, we added a comparison of the mean absolute error (MAE) and the root mean square error (RMSE) in the supplementary material. Table I shows that our proposed IHQS_p-CG algorithm outperforms other traditional algorithms, even surpassing the deep learning algorithms PDNet and FBPNet on both AAPM and Covid-19 datasets. It shows greater robustness to noise and missing range than the HQS-CG algorithm, with an average improvement of more than 0.5 dB for 180-degree and 120-degree observations, and over 1 dB for 90-degree and 60-degree observations.

Our proposed IHQS_p-Net outperforms other deep learning methods, especially in handling mixed noises. Our lightweight IHQS_p-Net* significantly improves upon the traditional IHQS_p-CG, indicating a lightweight CG initializer can enhance performance and convergence speed. Compared to FBPNet, IHQS_p-Net shows an average improvement of 2 dB, highlighting the benefits of integrating traditional algorithms with theoretical guarantees into deep networks. Meanwhile, Table I gives the running time for each algorithm to reconstruct a 512$\times$512 observation of 90 degrees.

Fig. 2 compares reconstructed images. FBP, BM3D, and HQS-CG show artifacts and missing structures due to lack of observation angle. Our IHQS_p-CG yields smoother results. Deep learning algorithms like PDNet and FBPNet produce smooth but over-smoothed results, blurring structural information. Both MetaInvNet and our IHQS_p-Net produce pleasing results, but IHQS_p-Net captures more details and avoids pseudo-connections seen in MetaInvNet. Notably, IHQS_p-Net* visually outperforms IHQS_p-CG, highlighting the importance of the CG initializer. The robustness analysis of all the algorithms is given in Fig. 3. By statistically analyzing the reconstructed images for all noise cases, observation angle cases and datasets, it can be found that the proposed IHQS_p-CG and IHQS_p-Net are robust to noise, observation angle, and data.

There is a key parameter $p$ in IHQS_p which controls the sparsity of the regular term. Different values of $p$ affect the reconstruction results of CT images [9]. Therefore, we need to discuss the different $p$-value effects. Fig. 4(a) shows the results of different $p$ values on two datasets, AAPM and Covid-19. From the figure, we can find that between 1 and 0.7, the reconstruction effect is enhanced as $p$ decreases. This reflects the advantage of the $L_{p}$ paradigm over $L_{1}$. However, lower$p$ is not better. As can be seen in the figure, for IHQS_p, the optimal $p=0.7$.

One of the contributions of IHQS_p is the introduction of inertial steps to accelerate the convergence of HQS. HQS_p-CG, in Table I, shows the numerical results without inertia step. It can be found that the introduction of the inertia step can effectively improve the algorithm’s performance, and the improvement becomes more obvious as the observation angle becomes sparser. Fig. 4(b) gives a comparison between IHQS_p and HQS_p-CG for iterative reconstruction when the 90-degree observation is corrupted by mixed Gaussian Poisson noise. PSNR and error $\frac{||u^{k}-u^{k-1}||^{2}}{||u^{k}||^{2}}$ are used as metrics. From Fig. 4(b), IHQS_p stops iterating at step 65, while HQS_p without inertial steps iterates until step 90. Meanwhile, HQS_p with more iterations does not outperform IHQS_p in PSNR. This reflects the importance of introducing inertial steps.

Table II gives the comparison of MAE and RMSE for all the compared algorithms. The results in Table II are similar to those in Table I. IHQ$S_{p}$-CG performs optimally in the traditional algorithm. IHQS_p-Net also performs best among the deep learning algorithms in aggregate. This shows the robust results of the proposed IHQS_p algorithm under different metrics.

Fig. 5 illustrates a qualitative comparison of the 60-degree observation resulting from the reconstruction. It can be noticed that MetaInvNet is prone to pseudo-generation. HQS-CG, on the other hand, is not able to complete the reconstruction details. Overall, IHQS_p-CG and IHQS_p-Net are good at preserving details in both traditional and deep learning algorithms.