2.1

Supervised Module

The supervised modules are trained sequentially. We set the loss function during training to be the root-mean-squared error (RMSE) to enforce alignment between the refined images and ground truth images. In the l-th parallel SUPER layer, the optimization problem for training the neural network is:

where G_θ(l) (·) denotes the neural network mapping in the l-th layer with parameters θ^(l), is the n-th input image from (l - 1)-th layer, is the corresponding regular-dose (reference) image or the training label. Note that the neural networks in different layers have different parameters.

2.2

Unsupervised Module

For the unsupervised module of each layer, we solve the following MBIR problem to reconstruct an image x ∈ ℝ^N_p from the corresponding noisy sinogram data y ∈ ℝ^N_d:

where W = diag{w_i} ∈ ℝ^{N_d ⨯ N_d} is a diagonal weighting matrix with the diagonal elements w_i being the estimated inverse variance of y_i, A ∈ ℝ ^{N_d ⨯ N_p} is the system matrix of the CT scan, L(Ax, y) is the data-fidelity term, penalty ℜ(x) is a (learning-based) regularizer, and the parameter β > 0 controls the noise and resolution trade-off.

In this work, we use the PWLS-ULTRA method to reconstruct an image x from noisy sinogram data y (measurements) with a union of pre-learned transforms . The image reconstruction is done through the following nonconvex optimization problem:

where is the reconstructed image by the unsupervised solver in the l-th layer, the operator P_j ∈ ℝ^{l⨯ N_p} extracts the j-th patch of l voxels of image x as P_jx, z_j is the corresponding sparse encoding of the image patch under a matched transform, and 𝒞_k, denotes the indices of patches grouped into the k-th cluster with transform Ω_k. Minimization over 𝒞_k indicates the computation of the cluster assignment of each patch. The regularizer ℜ includes an encoding error term and an ℓ₀ sparsity penalty counting the number of non-zero entries with weight ^γ2. The sparse encoding and clustering are computed simultaneously. We apply the alternating minimization method from⁵ (with inner iterations for updating x) on the above optimization problem. The algorithm also uses a different (potentially better) initialization in each parallel SUPER layer, which may benefit solving the involved nonconvex optimization problem.

2.3

Parallel SUPER Model

The main idea of the Parallel SUPER framework is to combine the supervised neural networks and iterative model-based reconstruction solvers in each layer. Define as an iterative MBIR solver with initial solution noisy sinogram data y and hyperparameter setting Γ to solve optimization problem (2). In the l-th layer, the parallel SUPER model is formulated as:

where λ is the nonnegative weight parameter for the neural network output in each layer and it is selected and fixed in all layers. Each parallel super layer can be thought of as a weighted average between a supervised denoised image and a reconstructed low-dose image from the unsupervised solver. Repeating multiple parallel SUPER layers simulates a fixed-point iteration to generate an ultimate reconstructed image.

The Parallel SUPER training algorithm based on (P0) is shown in Algorithm 1. The Parallel SUPER reconstruction algorithm is the same except that it uses the learned network weights in each layer.

3.1

Experiment Setup

In our experiments, we use the Mayo Clinics dataset established for “the 2016 NIH-AAPM-Mayo Clinic Low Dose CT Grand Challenge”.¹⁰ We choose 520 images from 6 of 10 patients in the dataset, among which 500 slices are used for training and 20 slices are used for validation. We randomly select 20 images from the remaining 4 patients for testing. We project the regular dose CT images x^⋆ to sinograms y by adding Poisson and additive Gaussian noise to them as follows:

where the original number of incident photons per ray is I₀ = 10⁴, the Gaussian noise variance is σ² = 25, and ϵ¹¹ is a small positive number to avoid negative measurement data when taking the logarithm.

We use the Michigan Image Reconstruction Toolbox to construct fan-beam CT geometry with 736 detectors ⨯ 1152 regularly spaced projection views, and a no-scatter mono-energetic source. The width of each detector column is 1.2858 mm, the source to detector distance is 1085.6 mm, and the source to rotation center distance is 595 mm. We reconstruct images of size 512 ⨯ 512 with the pixel size being 0.69 mm ⨯ 0.69 mm.

3.2

Parameter Settings

In the parallel SUPER model, we use FBPConvNet as the supervised module and PWLS-ULTRA as the unsupervised module. It takes about 10 hours for training the model for 10 layers in a GTX Titan GPU graphics processor. We train models for different values of the parameter λ (to then select an optimal value), including 0.1, 0.3, 0.5, 0.7, and 0.9. During the training of the supervised method, we ran 4 epochs (kept small to reduce overfitting risks) of the stochastic gradient descent (SGD) optimizer for the FBPConvNet module in each parallel SUPER layer. The training hyperparameters of FBPConvNet are set as follows: the learning rate decreases logarithmically from 0.001 to 0.0001; the batchsize is 1; and the momentum parameter is 0.99. The filters are initialized in the various networks during training with i.i.d. random Gaussian entries with zero mean and variance 0.005. For the unsupervised module, we have trained a union of 5 sparsifying transforms using 12 slices of regular-dose CT images (which are included in the 500 training slices). Then, we use the pre-learned union of 5 sparsifying transforms to reconstruct images with 5 outer iterations and 5 inner iterations of PWLS-ULTRA. In the training and reconstruction with ULTRA, we set the parameters β = 5 ⨯ 10³ and γ = 20. PWLS-EP reconstruction is used as the initialization of the input of network in the first layer.

We compare the proposed parallel SUPER model with the unsupervised method (PWLS-EP), standalone supervised module (FBPConvNet), standalone unsupervised module (PWLS-ULTRA), and the serial SUPER model. PWLS-EP is a penalized weighted-least squares reconstruction method with edge-preserving hyperbola regularization. For the unsupervised method (PWLS-EP), we set the parameters δ = 20 and β = 2¹⁵ and run 100 iterations to obtain convergent results. In the training of the standalone supervised module (FBPConvNet), we run 100 epochs of training to sufficiently learn the image features with low overfitting risks. In the standalone unsupervised module (PWLS-ULTRA), we use the pre-learned union of 5 sparsifying transforms to reconstruct images. We set the parameters β = 10⁴ and γ = 25, and run 1000 alternations with 5 inner iterations to ensure good performance. In the serial SUPER model, we run 4 epochs of training when learning the supervised modules (FBPConvNet), and we use the pre-learned union of 5 sparsifying transforms and set the parameters β = 5 ⨯ 10³, γ = 20 and μ = 5 ⨯ 10⁵ to reconstruct images with 20 alternations and 5 inner iterations for the unsupervised module (PWLS-ULTRA).

3.3

Results

To compare the performance quantitatively, we compute the RMSE in Hounsfield units (HU) and structural similarity index measure (SSIM)¹² for the reconstructed images. For a reconstructed image , RMSE is defined as, , where denotes the reference regular-dose image intensity at the j-th pixel location and N_p is the number of pixels.

We train the parallel SUPER framework with different choices of the parameter λ including 0.1, 0.3, 0.5, 0.7 and 0.9 to obtain the best choice. Fig. 2 shows the evolution of RMSE over layers for 20 validation slices with different λ choices. We can see that we obtain the best RMSE when λ = 0.3.

Figure 2:

RMSE (over 20 test slices) comparison with different choices of parameters.

We have conducted experiments on 20 test slices (slice 20, slice 50, slice 100, slice 150 and slice 200 of patient L067, L143, L192, L310) of the Mayo Clinic data. Table 1 shows the averaged image quality of 20 test images with different methods. From Table 1, we observe that Parallel SUPER significantly improves the image quality compared with the standalone methods. It also achieves 1.8 HU better average RMSE compared with Serial SUPER while its SSIM is comparable with Serial SUPER. Fig. 3 shows the reconstructions of L067 (slice 50) and L310 (slice 150) using PWLS-ULTRA, FBPConvNet, serial SUPER (FBPConvNet + PWLS-ULTRA), and parallel SUPER (FBPConvNet + PWLS-ULTRA), along with the references (ground truth). The Parallel SUPER scheme achieved the lowest RMSE and the zoom-in areas show that Parallel SUPER can reconstruct image details better.

Figure 3:

Reconstruction of slice 50 from patient L067 and reconstruction of slice 150 from patient L310 using various methods. The display window is [800, 1200] HU.

Table 1:

The mean RMSE and SSIM values for 20 test images with PWLS-EP, PWLS-ULTRA, FBPConvNet, Serial SUPER, and the proposed Parallel SUPER.