Neural networks-based regularization for large-scale medical image reconstruction

Clearly, the regularization term $\mathcal{R}({\mathbf{x}})$ significantly affects the quality and the characteristics of the solution ${\mathbf{x}}$. Here, we propose a generalized approach for solving high-dimensional inverse problems by the following three steps: First, an initial guess of the solution is provided by a direct reconstruction from the measured data, i.e. ${\mathbf{x}}_{\mathrm{ini}} = {\mathbf A}^{\dagger}{\mathbf{y}}$, where ${\mathbf{A}}^{\dagger}\colon Y \to X$ denotes some reconstruction operator. Then, a CNN is used to remove the noise or the artefacts from the direct reconstruction ${\mathbf{x}}_{\mathrm{ini}}$ in order to obtain another intermediate reconstruction ${\mathbf{x}_{\mathrm{CNN}}}$ which is used as a CNN-prior in a generalized Tikhonov functional

$$\begin{equation} F_{{\mathbf{y}},{\mathbf{x}_{\mathrm{CNN}}},\lambda} ({\mathbf{x}}) : = D({\mathbf{A}} {\mathbf{x}},{\mathbf{y}}) + \lambda \| {\mathbf{x}} - {\mathbf{x}_{\mathrm{CNN}}}\|_2^2\rightarrow \mathrm{\textrm{min}}. \end{equation} \tag{ 3 }$$

As a third and final step, the CNN-Tikhonov functional (3) is minimized resulting in the proposed CNN-based reconstruction.

Note that the regularization of the problem, i.e. obtaining the CNN-prior, is decoupled from the step of ensuring data-consistency of the solution via minimization of (3). This means that, in our case, the regularization term $\mathcal{R}({\mathbf{x}})$ in (2) depends on the data ${\mathbf{y}}$ and is of the form

$$\begin{equation} \mathcal{R}({\mathbf{x}}) = \mathcal{R}({\mathbf{x}},f_{\theta}({\mathbf{A}}^{\dagger} {\mathbf{y}})), \end{equation} \tag{ 4 }$$

where $f_{\Theta}$ denotes a function which decomposes the initial reconstruction ${\mathbf A}^{\dagger} {\mathbf{y}}$ into patches, processes them with a pre-trained CNN $u_{\Theta}$ with trainable parameters Θ and recomposes the patches in order to provide the CNN image-prior ${\mathbf{x}_{\mathrm{CNN}}}$. Therefore, for minimizing (3) we only have to apply $f_{\Theta}$ one single time. Thus, the training of the network $u_{\Theta}$ is facilitated since it only has to learn to map the manifold given by the initial reconstructions to the manifold of the ground truth images.

This also allows to use deeper and more sophisticated CNNs as the ones typically used in iterative networks. Given the high-dimensionality of the considered problems, network training is further carried out on sub-portions of the image samples, i.e. on patches or slices which are previously extracted from the images or volumes. This is motivated by the fact that in most medical imaging applications, one has typically access to datasets with only a relatively small number of subjects. The images or volumes of these subjects, on the other hand, are elements of a high-dimensional space. Therefore, one is concerned with the problem of having topologically sparse training data with only very few data points in the original high-dimensional image space. Working with sub-portions of the image samples increases the number of available data points and at the same time decreases its ambient dimensionality.

Suppose we have access to a finite set of N ground truth samples $({\mathbf{x}}_{k})_{k = 1}^N$ and corresponding initial estimates $({\mathbf{x}}_{\mathrm{ini},k})_{k = 1}^N$. We are in particular interested in the case where N is relatively small and the considered samples ${\mathbf{x}}_{k}$ have a relatively large size, which is the case for most medical imaging applications. For any sample ${\mathbf{x}} \in X$ we consider its decomposition in $N_{\mathbf{p},\mathbf{s}}$ possibly overlapping patches

$$\begin{equation} {\mathbf{x}} = \mathbf{W}_{\mathbf{p},\mathbf{s}}\ \sum_{j = 1}^{N_{\mathbf{p},\mathbf{s}}} ({\mathbf{R}}_j^{\mathbf{p},\mathbf{s}})^{{{\scriptstyle \boldsymbol{\mathsf{T}}}}} \, {\mathbf{R}}_j^{\mathbf{p},\mathbf{s}} \, {\mathbf{x}}, \end{equation} \tag{ 5 }$$

where ${\mathbf{R}}_j^{\mathbf{p},\mathbf{s}}$ and $({\mathbf{R}}_j^{\mathbf{p},\mathbf{s}})^{{{\scriptstyle \boldsymbol{\mathsf{T}}}}}$ extract and reposition the patches at the original position, respectively, and the diagonal operator $\mathbf{W}_{\mathbf{p},\mathbf{s}}$ accounts for weighting of regions containing overlaps. The entries of the tuples $\mathbf{p}$ and $\mathbf{s}$ specify the size of the patches and the strides in each dimension and therefore the number of patches $N_{\mathbf{p},\mathbf{s}}$ which are extracted from a single image.

We aim for improved estimates ${\mathbf{x}}_{\mathrm{CNN},k} = f_\theta( {\mathbf{x}}_{\mathrm{ini},k}) \approx {\mathbf{x}}_{k}$ via a trained network function $f_\theta$ to be constructed. Since the operator norm of $\mathbf{W}_{\mathbf{p},\mathbf{s}}$ is less or equal to one, by the triangle inequality, we can estimate the average error

$$\begin{equation} e_N: = \sum_{k = 1}^N \|{\mathbf{x}}_{k} - {\mathbf{x}}_{\mathrm{CNN},k} \|_2 \leq \sum_{k = 1}^N \sum_{j = 1}^{N_{\mathbf{p},\mathbf{s}}} \bigl\| {\mathbf{R}}_j^{\mathbf{p},\mathbf{s}} {\mathbf{x}}_{k} - {\mathbf{R}}_j^{\mathbf{p},\mathbf{s}} {\mathbf{x}}_{\mathrm{CNN},k} \bigr\|_2 = : e_{N,N_{\mathbf{p},\mathbf{s}}}. \end{equation} \tag{ 6 }$$

Inequality (6) suggests that it is beneficial estimating each patch of the sample ${\mathbf{x}}_{k}$ by a neural network $u_{\theta}$ applied to ${\mathbf{R}}_j^{\mathbf{p},\mathbf{s}} {\mathbf{x}}_{\mathrm{ini},k}$ rather than estimating the whole sample at once. The neural network $u_{\theta}$ is trained on a subset of pairs

$$\begin{eqnarray} &&\mathcal{D} = \Big\{\Big({\mathbf{R}}_j^{\mathbf{p},\mathbf{s}} ({\mathbf{x}}_{\mathrm{ini},k}), {\mathbf{R}}_j^{\mathbf{p},\mathbf{s}} ({\mathbf{x}}_{k} ) \Big) : (k,j) \in \mathcal{I}_{N,N_{\mathbf{p},\mathbf{s}}} \Big\}, \end{eqnarray} \tag{ 7 }$$

of all possible patches extracted from the N samples in the dataset, where $\mathcal{I}_{N,N_{\mathbf{p},\mathbf{s}}}: = \{1,\ldots,N\} \times \{1,\ldots,N_{\mathbf{p},\mathbf{s}} \}$. During training, we optimize the set of parameters θ to minimize the L₂-error between the estimated output of the patches and the corresponding ground truth patch by minimizing

$$\begin{equation} \mathcal{L}(\theta) = \frac{1}{N_{\mathrm{train}}} \sum_{(\mathbf{z}_{\mathrm{ini}}, \mathbf{z}) \in \mathcal{D} }\| u_{\theta}(\mathbf{z}_{\mathrm{ini}}) - \mathbf{z}\|_2^2 \,, \end{equation} \tag{ 8 }$$

where $N_{\mathrm{train}}$ is the number of training patches.

Denote by $f_{\theta}$ the composite function which decomposes a sample image or volume ${\mathbf{x}}$ into patches, applies a neural network $u_{\theta}$ to each patch, and reassembles the sample from them. This results in the proposed CNN-prior ${\mathbf{x}_{\mathrm{CNN}}}$ given by

$$\begin{equation} {\mathbf{x}_{\mathrm{CNN}}} : = f_{\theta}({\mathbf{x}}_{\mathrm{ini}}) = \mathbf{W}_{\mathbf{p},\mathbf{s}}\sum_j ({\mathbf{R}}_j^{\mathbf{p},\mathbf{s}})^{{{\scriptstyle \boldsymbol{\mathsf{T}}}}} ( u_{\theta}({\mathbf{R}}_j^{\mathbf{p},\mathbf{s}}({\mathbf{x}}_{\mathrm{ini}}) ) ), \end{equation} \tag{ 9 }$$

where ${\mathbf{x}}_{\mathrm{ini}} = {\mathbf{A}}^\dagger {\mathbf{y}}$ is the initial reconstruction obtained from the measured data.

Remark. Since in (6), the right-hand-side is an upper bound of the left-hand-side, inequality (6) guarantees that the set of parameters found by minimizing (8) is also suitable for obtaining the prior ${\mathbf{x}_{\mathrm{CNN}}}$. Therefore, $u_{\theta}$ is powerful enough to deliver a CNN-prior to regularize the solution of (3). Figure 1 illustrates the process of extracting patches from a volume using the operator ${\mathbf{R}}_j^{\mathbf{p},\mathbf{s}}$, processing it with a neural network $u_{\theta}$ and repositioning it at the original position using the transposed operator $({\mathbf{R}}_j^{\mathbf{p},\mathbf{s}})^{{{\scriptstyle \boldsymbol{\mathsf{T}}}}}$. The example is shown for a 2D cine MR image sequence.

Zoom In Zoom Out Reset image size

After having found the CNN prior (9), as a final reconstruction step, the optimality condition for minimization problem (3) is solved with an iterative method dependent on the specific application. Minimizing this functional increases data-consistency compared to ${\mathbf{x}_{\mathrm{CNN}}}$ since the discrepancy to the measured data is used again. The solution of (3) is then the final CNN-based reconstruction of the proposed method. Algorithm 1 summarizes the proposed three-step reconstruction scheme.

Note that the regularizer $\mathcal{R}({\mathbf{x}}) = \| {\mathbf{x}} - {\mathbf{x}_{\mathrm{CNN}}} \|_2^2$ is strongly convex. Therefore, if the discrepancy term $D({\mathbf{A}} {\mathbf{x}}, {\mathbf{y}})$ is convex, then the Tikhonov functional (3) is strongly convex and can be efficiently minimized by most gradient based iterative schemes including Landweber's iteration and Conjugate Gradient type methods. The specific strategy for minimizing (3) depends on the considered application. In the case of an ill-conditioned inverse problem with noisy measurements, it might be beneficial that (3) is only approximately minimized. For example, for the case of low-dose CT, early stopping of the Landweber iteration is applied as additional regularization method due to the semi-convergence property of the Landweber iteration (Strand 1974). Such a strategy is used for the numerical results presented in section 4.

Another benefit of our approach is that minimization of the Tikhonov functional (3) corresponds to convex variational regularization with a quadratic regularizer. Therefore, one can use well known stability, convergence and convergence rates results (Engl et al 1996, Scherzer et al 2009, Grasmair 2010, Li et al 2020). Consequently, opposed to most existing neural network based reconstruction algorithms, the proposed framework is build on the solid theoretical fundament for regularizing inverse problems. As an example of such results we have the following theorem. Convergence of CNN-based regularization

Theorem 1. Let ${\mathbf{A}} \colon X \to Y$ be linear, ${\mathbf{x}_{\mathrm{CNN}}} \in X$, ${\mathbf{y}}_0 \in {\mathbf{A}}(X)$, and ${\mathbf{y}}_\delta \in Y$ satisfy $\| {\mathbf{y}}_\delta -{\mathbf{y}}_0 \| \leq \delta$. Then the following hold:

• 1.  
  For all δ, λ > 0, the quadratic Tikhonov functional

  $$\begin{equation} F_{{\mathbf{y}}_\delta, {\mathbf{x}_{\mathrm{CNN}}},\lambda} ({\mathbf{x}}) : = \| {\mathbf{A}} {\mathbf{x}}-{\mathbf{y}}_\delta \|_2^2 + \lambda \| {\mathbf{x}} - {\mathbf{x}_{\mathrm{CNN}}}\|_2^2 \end{equation} \tag{ 10 }$$

  has a unique minimizer ${\mathbf{x}}_{\delta,\lambda}$.
• 2.  
  The exact data equation ${\mathbf{A}} {\mathbf{x}} = {\mathbf{y}}_0$ has a unique ${\mathbf{x}_{\mathrm{CNN}}}$-minimizing solution ${\mathbf{x}}_0 \in \arg\min \{\| {\mathbf{x}} - {\mathbf{x}_{\mathrm{CNN}}}\|^2 : {\mathbf{A}} {\mathbf{x}} = {\mathbf{y}}_0 \}$.
• 3.  
  If the parameter choice λ = λ(δ) satisfies λ, δ²/λ → 0 as δ → 0, then $\lim_{\delta \to 0} \| {\mathbf{x}}_0 - {\mathbf{x}}_{\delta,\lambda} \| = 0$.

Proof. The change of variables

• $\bar {\mathbf{x}} : = {\mathbf{x}} - {\mathbf{x}_{\mathrm{CNN}}}$
• $\bar {\mathbf{y}}_0 : = {\mathbf{y}}_0 - {\mathbf{A}} {\mathbf{x}_{\mathrm{CNN}}}$
• $\bar {\mathbf{y}}_\delta : = {\mathbf{y}}_\delta - {\mathbf{A}} {\mathbf{x}_{\mathrm{CNN}}}$

reduces (10) to standard Tikhonov regularization $\| {\mathbf{A}} \bar {\mathbf{x}} - \bar {\mathbf{y}}_\delta \|^2 + \lambda \| \bar {\mathbf{x}} \|_2^2 \to \min_{\mathbf{x}}$ for the inverse problem ${\mathbf{A}} \bar {\mathbf{x}} = \bar {\mathbf{y}}_0$. Therefore, Items (1.) – (3.) follow from standard results that can be found for example in (Engl et al 1996, section 5).

Algorithm 1: Proposed CNNs-based large-scale image reconstruction algorithm.
Data: trained network $u_{\theta}$, function $f_{\theta}$, noisy or incomplete measured data ${\mathbf{y}}$, regularization parameter λ > 0
Output: reconstruction ${\mathbf{x}}_{\mathrm{REC}}$
1) ${\mathbf{x}}_{\mathrm{ini}} \gets {\mathbf{A}}^{\dagger} {\mathbf{y}}$
2) ${\mathbf{x}_{\mathrm{CNN}}} \gets f_{\theta}({\mathbf{x}}_{\mathrm{ini}})$
3) ${\mathbf{x}}_{\mathrm{REC}} \gets \arg\min_{{\mathbf{x}}} D({\mathbf{A}} {\mathbf{x}}, {\mathbf{y}}) + \lambda \| {\mathbf{x}} - {\mathbf{x}_{\mathrm{CNN}}} \|_2^2$
Return ${\mathbf{x}}_{\mathrm{REC}}$

Theorem 1 also holds in the infinite-dimensional setting (Engl et al 1996, section 5) reflecting the stability of the proposed CNN regularization. Similar results hold for non-linear problems and general discrepancy measures (Grasmair 2010). Moreover, one can derive quantitative error estimates similar to Scherzer et al (2009), Grasmair (2010), Li et al (2020). Such theoretical investigations, however, are beyond the scope of this paper.

Here we applied our method to image reconstruction in undersampled 2D radial cine MRI. Typically, MRI is performed using multiple receiver coils and therefore, the inverse problem is given by

$$\begin{equation} {\mathbf{E}_I} {\mathbf{x}} = {\mathbf{y}_I}, \end{equation} \tag{ 11 }$$

where ${\mathbf{x}} \in \mathbb{C}^N$ with $N = N_x\cdot N_y\cdot N_t$ is an unknown complex-valued image sequence. The encoding operator ${\mathbf{E}_I}$ is given by ${\mathbf{E}_I} = \mathbf{S}\circ {\mathbf E} \circ {\mathbf C}$ where

$$\begin{eqnarray} {\mathbf C} & = [ {\mathbf C}_1,\ldots, {\mathbf C}_{n_c}]^{{{\scriptstyle \boldsymbol{\mathsf{T}}}}}, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} {\mathbf E} & = \mathrm{diag}( {\mathbf F},\ldots, {\mathbf F} ), \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} \mathbf{S} & = \mathrm{diag}( {\mathbf{S}_I},\ldots, {\mathbf{S}_I}). \end{eqnarray} \tag{ 14 }$$

Here, ${\mathbf C}_i$ denotes the ith coil sensitivity map, n_c is the number of coil-sensitivity maps, ${\mathbf F}$ the 2D frame-wise operator and ${\mathbf{S}_I}$ with $I \subset J = \{1,\ldots,N_{\mathrm{rad}}\}$, $|I|: = m \leq N_{\mathrm{rad}}$, a binary mask which models the undersampling process of the $N_{\mathrm{rad}}$ Fourier coefficients sampled on a radial grid. The vector ${\mathbf{y}_I} \in \mathbb{C}^M$ with M = m · n_c corresponds to the measured data. Here, we sampled the k-space data along radial trajectories chosen according to the golden-angle method (Winkelmann et al 2006). Note that problem (11) is mainly ill-posed not due to the presence of noise in the acquisition, but because the data acquisition is accelerated and hence only a fraction of the required measurements is acquired.

If we assume a radial data-acquisition grid, problem (11) is a large-scale inverse problem mainly because of two reasons. First, the measurement vector ${\mathbf{y}_I}$ corresponds to n_c copies of the Fourier encoded image data multiplied by the corresponding coil sensitivity map. Second, the adjoint operator ${\mathbf{E}_I}^{{\scriptstyle \boldsymbol{\mathsf{H}}}}$ consists of two computationally demanding steps. The radially acquired k-space data is first properly re-weighted and interpolated to a Cartesian grid, for example by using Kaiser-Bessel functions (Rasche et al 1999). Then, a 2D inverse Fourier operation is applied to the image of each cardiac phase and the final image sequence is obtained by weighting the images from each estimated coil-sensitivity map and combining them to a single image sequence. We refer to the reconstruction obtained by ${\mathbf{x}_I} = {\mathbf{E}_I}^{{\scriptstyle \boldsymbol{\mathsf{H}}}}{\mathbf{y}_I}$ as the non-uniform fast Fourier-transform (NUFFT) reconstruction. Therefore, in radial multi-coil MRI, the measured k-space data is high-dimensional and the application of the encoding operators ${\mathbf{E}_I}$ and ${\mathbf{E}_I}^{{\scriptstyle \boldsymbol{\mathsf{H}}}}$ is further more computationally demanding than sampling on a Cartesian grid, see e.g Smith et al (2019). This makes the construction of cascaded networks which also process the k-space data (Han et al 2019) or repeatedly employ the forward and adjoint operators (Schlemper et al 2018, Qin et al 2018) computationally challenging. Therefore, the separation of the regularization given by the CNNs from the data-consistency step is necessary in this case due to computational reasons.

As proposed in section 2, we solve a regularized version of problem (11) by minimizing

$$\begin{equation} F_{{\mathbf{y}_I}, {\mathbf{x}_{\mathrm{CNN}}}, \lambda} ({\mathbf{x}}) = \| {\mathbf{E}_I} {\mathbf{x}} - {\mathbf{y}_I} \|_2^2 + \lambda \| {\mathbf{x}} - {\mathbf{x}_{\mathrm{CNN}}} \|_2^2, \end{equation} \tag{ 15 }$$

where ${\mathbf{x}_{\mathrm{CNN}}}$ is obtained a-priori by using an already trained network. For this example, for obtaining the CNN-prior ${\mathbf{x}_{\mathrm{CNN}}}$, we adopted the XT,YT approach presented in Kofler et al (2019) which has been reported to achieve comparable results to the 3D U-net (Hauptmann et al 2019). Using the XT,YT approach has the main advantage of only using 2D convolutional layers while still exploiting the spatio-temporal information of the cine MR images. The XT,YT approach works by first extracting all x, t- and y, t-slices from the image, subsequently processing them with a modified version of the 2D U-net and then properly reassembling them to obtain an estimate of the artefact-free image. Since the network only has to learn the manifold of the x, t- and y, t-spatio-temporal slices, the approach was shown to be applicable even when only little training data is available. Since the XT,YT method was previously introduced to only process real-valued data (i.e. the magnitude images), we followed a similar strategy by processing the real and imaginary parts of the image sequences separately but using the same real-valued network $u_{\theta}$. This further increases the amount of training data by a factor of two. More precisely, let ${\mathbf{R}}^{xt}_j$ and ${\mathbf{R}}^{yt}_j$ denote the operators which extract the jth two-dimensional spatio-temporal slices in xt- and yt-direction from a 3D volume $({\mathbf{R}}^{xt}_j)^{{\scriptstyle \boldsymbol{\mathsf{T}}}}$ and $({\mathbf{R}}^{yt}_j)^{{\scriptstyle \boldsymbol{\mathsf{T}}}}$ their respective transposed operations which reposition the spatio-temporal slices at their original position.

By $u_{\theta}$ we denote a 2D U-net as the one described in Kofler et al (2019) which is trained on spatio-temporal slices, i.e. on a dataset of pairs which consist of the spatio-temporal slices in xt- and yt-direction of both the real and imaginary parts of the complex-valued images. The network $u_{\theta}$ was trained to minimize the L₂-error between the ground truth image and the estimated output of the CNN. Our dataset consists of radially acquired 2D cine MR images from n = 19 subjects (15 healthy volunteers and 4 patients with known cardiac dysfunction) with 30 images covering the cardiac cycle. The ground truth images were obtained by kt-SENSE reconstruction using $N_{\theta} = 3400$ radial lines. We retrospectively generated the radial k-space data ${\mathbf{y}_I}$ by sampling the k-space data along $N_{\theta} = 1130$ radial spokes using n_c = 12 coils. Note that sampling $N_{\theta} = 3400$ already corresponds to an acceleration factor of approximately ∼ 3 and therefore, $N_{\theta} = 1130$ corresponds to an accelerated data-acquisition by an approximate factor of ∼ 9. The forward and the adjoint operators ${\mathbf E}_I$ and ${\mathbf E}_I^{{\scriptstyle \boldsymbol{\mathsf{H}}}}$ were implemented using the ODL library (Adler et al 2017). The complex-valued CNN-regularized image sequence ${\mathbf{x}_{\mathrm{CNN}}}$ was obtained by

$$\begin{align} {\mathbf{x}_{\mathrm{CNN}}} & = f_{\theta}({\mathbf{x}_I}) \nonumber\\ \nonumber & = \frac{1}{2} \Big[ \sum_j ({\mathbf{R}}^{xt}_j)^{{{\scriptstyle \boldsymbol{\mathsf{T}}}}} \big(u_{\theta}( {\mathbf{R}}^{xt}_j (\mathrm{\textrm{Re}}\,{\mathbf{x}_I})) \big) \nonumber + ({\mathbf{R}}^{yt}_j)^{{{\scriptstyle \boldsymbol{\mathsf{T}}}}} \big(u_{\theta}( {\mathbf{R}}^{yt}_j (\mathrm{\textrm{Re}}\,{\mathbf{x}_I})) \big) \\ & \quad + \mathrm{i}\, \Big(({\mathbf{R}}^{xt}_j)^{{{\scriptstyle \boldsymbol{\mathsf{T}}}}} \big(u_{\theta}( {\mathbf{R}}^{xt}_j (\mathrm{\textrm{Im}}\, {\mathbf{x}_I}))\big)\Big) \nonumber + \mathrm{i}\, \Big(({\mathbf{R}}^{yt}_j)^{{{\scriptstyle \boldsymbol{\mathsf{T}}}}} \big(u_{\theta}( {\mathbf{R}}^{yt}_j (\mathrm{\textrm{Im}}\, {\mathbf{x}_I}))\big)\Big)\Big ]. \end{align} \tag{ 16 }$$

Given ${\mathbf{x}_{\mathrm{CNN}}}$, functional (15) was minimized by setting its derivative with respect to ${\mathbf{x}}$ to zero and applying the pre-conditioned conjugate gradient (PCG) method to iteratively solve the resulting system. PCG was used to solve the system ${\mathbf H} {\mathbf{x}} = \mathbf{b}$ with

$$\begin{eqnarray} &&{\mathbf H} = {\mathbf{E}_I}^{{\scriptstyle \boldsymbol{\mathsf{H}}}}{\mathbf{E}_I} + \lambda\, {\mathbf I}, \nonumber \\ &&\mathbf{b} = {\mathbf{x}_I} + \lambda\, {\mathbf{x}_{\mathrm{CNN}}}. \end{eqnarray} \tag{ 17 }$$

Since the XT,YT method gives access to a large number of training samples, training the network $u_{\Theta}$ for 12 epochs was sufficient. The CNN was trained by minimizing the L₂-norm of the error between labels and output by using the Adam optimizer (Kingma and Ba 2014). We split our dataset in 12/3/4 subjects for training, validation and testing and performed a 4-fold cross-validation. For the experiment, we performed $n_{\mathrm{iter}} = 16$ subsequent iterations of PCG and empirically set λ = 0.1. Note that due to strong convexity, (15) has a unique minimizer and solving system (17) yields the desired minimizer. The obtained results can be found in Subsection 4.1.

The current generation of CT scanners performs the data-acquisition by emitting x-rays along trajectories in the form of a cone-beam for each angular position of the scanner. Therefore, for each angle φ of the rotation, one obtains an x-ray image which is measured by the detector array and thus, the complete sinogram data can be identified with a 3D array of shape $(N_{\phi},N_{r_x},N_{r_y})$. Thereby, $N_{\phi}$ corresponds to the number of angles the rotation of the scanner is discretized by and $N_{r_x}$ and $N_{r_y}$ denote the number of elements of the detector array. The values of these parameters vary from scanner to scanner but are in the order of $N_{\phi}\approx 1000$ for a full rotation of the scanner and $N_{r_x} \times N_{r_y} \approx 320 \times 800$ for a 320-row detector array, which is for example used for cardiac CT scans (Dewey et al 2009). The volumes obtained from the reconstructions are typically given by an in-plane number of pixels of $N_x \times N_y = 512 \times 512$ and varying number of slices N_z, dependent on the specific application. For this example, we consider a similar set-up as in Adler and Öktem (2017). The non-linear problem is given by

$$\begin{equation} {\mathbf{y}_{\eta}} = \mathbf{T} {\mathbf{x}} + {\boldsymbol{\eta}} = p \, \mathrm{\textrm{exp}} \{- \mu\, {\mathbf{R}} {\mathbf{x}}\} + {\boldsymbol{\eta}}, \end{equation} \tag{ 18 }$$

where p denotes the average number of photons per pixel, µ is the linear attenuation coefficient of water, ${\mathbf{R}}$ corresponds to the discretized version of a ray-transform with cone-beam geometry and the vector ${\boldsymbol{\eta}}$ denotes the Poisson-distributed noise in the measurements. Following our approach, we are interested in solving

$$\begin{equation} F_{{\mathbf{y}_{\eta}}, {\mathbf{x}_{\mathrm{CNN}}}, \lambda} ({\mathbf{x}}) = D_{\mathrm{KL}}(\mathbf{T}{\mathbf{x}}, {\mathbf{y}_{\eta}}) + \lambda \| {\mathbf{x}} - {\mathbf{x}_{\mathrm{CNN}}} \|_2^2 \rightarrow \mathrm{\textrm{min}}, \end{equation} \tag{ 19 }$$

where $D_{\mathrm{KL}}$ denotes the Kullback-Leibler divergence which corresponds to the log-likelihood function for Poisson-distributed noise. According to the previously introduced notation, the prior ${\mathbf{x}_{\mathrm{CNN}}}$ is given by ${\mathbf{x}_{\mathrm{CNN}}} = f_{\theta}({\mathbf{x}_{\eta}})$, where $f_{\theta}$ denotes a CNN-based processing method with trainable parameters θ and ${\mathbf{x}_{\eta}} = {\mathbf{R}}^{\dagger}(-\mu^{-1}\mathrm{\textrm{ln}}(p^{-1}{\mathbf{y}_{\eta}}))$ with ${\mathbf{R}}^{\dagger}$ being the filtered back-projection (FBP) reconstruction.

Since our object of interest ${\mathbf{x}}$ is a volume, it is intuitive to choose a NN which involves 3D convolutions in order to learn the filters by exploiting the spatial correlation of adjacent voxels in x-, y- and z-direction. In this particular case, $u_{\theta}$ denotes a 3D U-net similar to the one presented in Hauptmann et al (2019). Due to the large dimensionality of the volumes ${\mathbf{x}}$, the network $u_{\theta}$ cannot be applied to the whole volume. Instead, following our approach, the volume was divided into patches to which the network $u_{\theta}$ is applied. Therefore, the output ${\mathbf{x}_{\mathrm{CNN}}}$ was obtained as described in (9), where $u_{\theta}$ operates on 3D patches given by the vector $\mathbf{p} = (128,128,16)$, which denotes the maximal size of 3D patches which we were able to process by a 3D U-net. The strides used for the extraction and the reassembling of the volumes used in (9) is empirically chosen to be $\mathbf{s} = (16,16,8)$.

Training of the network $u_{\theta}$ was performed on a dataset of pairs according to (7), where we retrospectively generated the measurements ${\mathbf{y}_{\eta}}$ by simulating a low-dose scan on the ground truth volumes. For the experiment, we used 16 CT volumes from the randomized DISCHARGE trial (Napp et al 2017) which we cropped to a fixed size of 512 × 512 × 128. The simulation of the low-dose scan was performed as described in Adler and Öktem (2017) by setting p = 10 000 and µ = 0.02. The operator ${\mathbf{R}}$ is assumed to perform $N_{\phi} = 1000$ projections which are measured by a detector array of shape $N_{r_x} \times N_{r_y} = 320 \times 800$. For the implementation of the operators, we used the ODL library (Adler et al 2017). The source-to-axis and source-to-detector distances were chosen according to the DICOM files. Since the dataset is relatively small, we performed a 7-fold cross-validation where for each fold we split the dataset in 12 patients for training, 2 for validation and 2 for testing. The number of training samples $N_{\mathrm{train}}$ results from the number of patches times the number of volumes contained in the training set. We trained the network $u_{\theta}$ for 115 epochs by minimizing the L₂-norm of the error between labels and outputs. For training, we used the Adam optimizer (Kingma and Ba 2014). With the described configuration of $\mathbf{p}$ and $\mathbf{s}$, the resulting number of patches to be processed in order to obtain the prior ${\mathbf{x}_{\mathrm{CNN}}}$ is therefore given by $N_{\mathbf{p},\mathbf{s}} = 9\,375$. In this example, the solution ${\mathbf{x}}_{\mathrm{REC}}$ to problem (19) was then obtained by performing $n_{\mathrm{iter}} = 4$ iterations of Landweber's method where we further used the filtered-back projection ${\mathbf{R}}^{\dagger}$ as a left-preconditioner to accelerate the convergence of the scheme. For the derivation of the gradient of (19) with respect to ${\mathbf{x}}$, we refer to Adler and Öktem (2017). The regularization parameter was empirically set to λ = 1. The results can be found in subsection 4.2.

Here we discuss the methods of comparison in more detail and report the times needed to process and reconstruct the images or volumes. The data-discrepancy term D(·, ·) was again chosen according to the considered examples as previously discussed. The TV-minimization approach used for comparison is given by solving

$$\begin{equation} D({\mathbf{A}} {\mathbf{x}} , {\mathbf{y}}) + \lambda \| \mathbf{G}{\mathbf{x}} \|_1 \to \min_{{\mathbf{x}}} \end{equation} \tag{ 20 }$$

where $\mathbf{G}$ denotes the discretized version of the isotropic first order finite differences filter in all three dimensions, i.e. in x-, y- and t-direction for the MR example and in x-, y- and z-direction for the CT example. The solution of problem (20) was obtained by introducing an auxiliary variable $\mathbf{z}$ and alternating between solving for ${\mathbf{x}}$ and $\mathbf{z}$. For the solution of one of the sub-problems, an iterative shrinkage method was used, see Chambolle (2005) for more details. The second resulting sub-problem was solved by iteratively solving a system of linear equations, either by Landweber for the CT example or by PCG for the MRI example, as mentioned before.

The dictionary learning-based method used for comparison is given by the solution of the problem

$$\begin{equation} D( {\mathbf{A}} {\mathbf{x}}, {\mathbf{y}}) + \lambda \| {\mathbf{x}} - {\mathbf{x}}_{\mathrm{DIC}}\|_2^2 \to \min_{{\mathbf{x}}} \end{equation} \tag{ 21 }$$

where, in contrast to our proposed method, ${\mathbf{x}}_{\mathrm{DIC}}$ was obtained by the patch-wise sparse approximation of the initial image estimate using an already trained dictionary $\mathbf{D}$. Therefore, using a similar notation as in (9), the prior ${\mathbf{x}}_{\mathrm{DIC}}$ is given by

$$\begin{equation} {\mathbf{x}}_{\mathrm{DIC}} = \mathbf{W}_{\mathbf{p},\mathbf{s}}\sum_j ({\mathbf{R}}_j^{\mathbf{p},\mathbf{s}})^{{{\scriptstyle \boldsymbol{\mathsf{T}}}}}\, \mathbf{D} \boldsymbol{\gamma}_j, \end{equation} \tag{ 22 }$$

where the dictionary $\mathbf{D}$ was previously trained by 15 iterations of the iterative thresholding and K residual means algorithm (ITKRM) (Schnass 2018) on a set of ground truth images which were given by the high-dose images for the CT example and the kt-SENSE reconstructions from $N_{\theta} = 3400$ radial lines for the MRI example. Here again, the dictionary $\mathbf{D}$ operates on three-dimensional patches, for the CT example/ x, y, z)-patches and the MR example (x, y, t)-patches. Note that for each fold, for training the dictionary $\mathbf{D}$, we only used the data which we included in the training set for our method. This means we trained a total of seven dictionaries for the CT example and four dictionaries for the MRI example. For each iteration of ITKRM, we randomly selected a subject to extract 10 000 3D training patches. The corresponding sparse codes $\boldsymbol{\gamma}_j$ were then obtained by solving

$$\begin{equation} \mathrm{\textrm{min}}_{\{\boldsymbol{\gamma}_j\}_j} \sum_j \| ({\mathbf{R}}_j^ {\mathbf{p},\mathbf{s}}) {\mathbf{x}}_{\mathrm{ini}} - \mathbf{D} \boldsymbol{\gamma}_j\|_2^2 + \| \boldsymbol{\gamma}_j\|_0, \end{equation} \tag{ 23 }$$

which is a sparse coding problem and was solved using orthogonal matching pursuit (OMP) (Tropp and Gilbert 2007). Thereby, the image ${\mathbf{x}}_{\mathrm{ini}}$ corresponds to either the FBP-reconstruction ${\mathbf{x}}_{\eta}$ for the CT example or to the NUFFT-reconstruction ${\mathbf{x}_I}$ for the MRI example. In both cases, we used patches of shape given by $\mathbf{p} = (4,4,4)$ and strides given by $\mathbf{s} = (2,2,2)$. The number of atoms K and the sparsity levels were set to K = 4 · d, with d = 4 · 4 · 4 and S = 16. Note that, in contrast to Xu et al (2012) and Wang et al (2004), Caballero et al (2014), the dictionary and the sparse codes were not learned during the reconstruction, as the sparse coding step of all patches would be too time consuming for very large-scale inverse problems, such as the CT example. Instead, the dictionary and the sparse codes were used to generate the prior ${\mathbf{x}}_{\mathrm{DIC}}$ which makes the method also more similar and comparable to ours. The parameter λ is set as previously stated in the manuscript, depending on the considered example.

For the evaluation of the reconstructions we report the normalized root mean squared error (NRMSE) and the peak signal-to-noise ratio (PSNR) as error-based measures and the structural similarity index measure (SSIM) (Wang et al 2004) and the Haar Wavelet-based perceptual similarity index measure (HPSI) (Reisenhofer et al 2018) as image-similarity-based measures. The reported statistics were obtained by calculating the measures of the images in the xy-plane and averaging them over the different folds.

Figure 2 shows an example of the results obtained with our proposed method. Figure 2(A) shows the initial NUFFT-reconstruction ${\mathbf{x}_I}$ obtained from the undersampled k-space data ${\mathbf{y}_I}$. The CNN-prior ${\mathbf{x}_{\mathrm{CNN}}}$ obtained by the XT,YT network can be seen in figure 2(B) and (shows) a strong reduction of undersampling artefacts but also blurring of small structures as indicated the yellow arrows. The CNN-prior ${\mathbf{x}_{\mathrm{CNN}}}$ is then used as a prior in functional (15) which is subsequently minimized in order to obtain the solution ${\mathbf{x}}_{\mathrm{REC}}$ which can be seen in figure 2(C). Figure 2(D) shows the kt-SENSE reconstruction from the complete sampling pattern using $N_{\theta} = 3400$ radial spokes for the acquisition. From the point-wise error images, we clearly see that the NRMSE is further reduced after performing the further iterations to minimize the CNN-prior-regularized functional. Further, fine details are recovered as can be seen from the yellow arrows in figure 2(C).

Zoom In Zoom Out Reset image size

Figure 3 shows a comparison of all different reported methods. As can be seen from the point-wise error in figure 3(B), the TV-minimization (Block et al 2007) method was able to eliminate some artefacts but less accurately compared to both learning-based methods, see figure 3(C) and (figure) 3(D). Table 1 lists the obtained quantitative measures for all methods averaged over the 4 different folds. From table 1, we see that the DIC method yielded better results than TV with respect to all reported measures. Our proposed solution ${\mathbf{x}}_{\mathrm{REC}}$ further surpassed the dictionary learning-based method, by additionally increasing the PSNR and SSIM by approximately 3 dB and 0.04, respectively. The difference with respect to HPSI, on the other hand, is relatively small. Our method also reduced the NRMSE by about 0.014 compared to the DIC method. In addition, from table 1, we see that for this example, even though processing the initial NUFFT-reconstruction with a CNN improved image quality with respect to all reported measures, further iterations to minimize the CNN-prior regularized functional increased data-consistency and additionally improved the PSNR, SSIM, HPSI and NRMSE. In fact, the statistics of the CNN-prior show that only post-processing the initial NUFFT-reconstruction leads to results which are inferior to the DIC method with respect to all reported measures.

Zoom In Zoom Out Reset image size

Table 1. Quantitative measures for the 2D radial cine MRI example. The measures are obtained as averages over the four different folds.

	NUFFT	${\mathbf{x}_{\mathrm{CNN}}}$	${\mathbf{x}}_{\mathrm{REC}}$	TV	DIC
PSNR	36.802 3	42.564 7	48.775 2	41.696 8	45.474 3
NRMSE	0.122 8	0.061 2	0.030 2	0.069 3	0.044 2
SSIM	0.664 9	0.787 6	0.952	0.863 5	0.917 5
HPSI	0.967 9	0.991 0	0.998 5	0.987 8	0.995 9

Figure 4 shows all the intermediate results obtained with the proposed method. Figure 4(A) shows the initial FBP-reconstruction which is contaminated by noise. The FBP-reconstruction was then processed using the function $f_{\theta}$ described in (9) to obtain the prior ${\mathbf{x}_{\mathrm{CNN}}}$ which can be seen in figure 4(B). From the point-wise error, we see that patch-wise post-processing with the 3D U-net removed a large portion of the noise resulting from the low-dose acquisition. Solving problem (19) increases data-consistency since we make use of the measured data ${\mathbf{y}_{\eta}}$. Note that in contrast to the previous example of undersampled radial MRI, the minimization of the functional increased data-consistency of the solution but also contaminated the solution with noise, since the measured data is noisy due to the simulated low-dose scan protocol. Table 2 summarizes the obtained quantitative measures for all intermediate reconstructions of our approach as well as for the TV and the DIC method. In the first three columns of table 2 we see the results obtained for all three intermediate reconstructions of our proposed scheme. The reconstruction metrics improved substantially from the FBP-reconstruction to the estimated prior ${\mathbf{x}_{\mathrm{CNN}}}$. The difference in terms of PSNR was almost 10 dB, while the NRMSE decreased by approximately 0.11. Further, the similarity measures SSIM and HPSI were increased by about 0.14 and 0.04, respectively. Finally, the estimated solution given by ${\mathbf{x}}_{\mathrm{REC}}$ which was obtained by performing $n_{\mathrm{iter}} = 4$ iterations of Landweber to minimize (19) showed a slight decrease in PSNR and NRMSE which is related to the use of the noisy-measured data. However, fine diagnostic details as the coronary arteries are still visible in the prior ${\mathbf{x}_{\mathrm{CNN}}}$ and in the solution ${\mathbf{x}}_{\mathrm{REC}}$ as indicated by the yellow arrows. SSIM slightly increased while HPSI stayed approximately the same.

Zoom In Zoom Out Reset image size

Table 2. Quantitative measures for the 3D low-dose CT example. The measures are obtained as averages over the seven different folds.

	FBP	${\mathbf{x}_{\mathrm{CNN}}}$	${\mathbf{x}}_{\mathrm{REC}}$	TV	DIC
PSNR	30.005 2	40.354 6	39.626 4	33.946	34.780 7
NRMSE	0.165 7	0.049 8	0.053 8	0.105 1	0.093 8
SSIM	0.425	0.575 5	0.581 3	0.498 5	0.546 5
HPSI	0.939 4	0.982 1	0.981 9	0.950 3	0.958 1

Figure 5 shows a comparison of images obtained by the different reconstruction methods. In figure 5(A), we see again the FBP-reconstruction obtained from the noisy data. Figure 5(B) shows the result obtained by the TV-minimization method which removed some of the noise as can be taken from the point-wise error image. The result obtained by the DIC method can be seen in figure 5(C) which further reduced image noise compared to the TV method and surpasses TV with respect to the reported statistics, as can be seein in table 2. Finally, figure 5(D) shows the solution ${\mathbf{x}}_{\mathrm{REC}}$ obtained with our proposed scheme and figure 5(E) shows the ground truth image. The reconstruction using the CNN output as a prior further increased the SSIM of the final result and is visually more similar to the ground truth image. The NRMSE and PSNR on the other hand were slightly reduced due to the use of the noisy measured data. HPSI remained approximately the same as can be seen from table 2.

Zoom In Zoom Out Reset image size

Table 3 summarizes the times for the different components of the reconstructions using all different approaches for both examples. The abbreviations 'SHRINK' and 'LS1' stand for 'shrinkage' and 'linear system - one iteration' and denote the times which are needed to apply the iterative shrinkage method for the TV approach and to solve the sub-problems which are solved using iterative schemes, respectively.

Table 3. Reconstruction and processing times for the different methods for one 3D CT volume and a 2D cine MR image sequence.

		3D Low Dose CT	2D Radial Cine MRI
${\mathbf{A}}^{\dagger} {\mathbf{y}_{\eta}}$		≈ 23 s (FBP)	≈ 11 s (NUFFT)
TV	$\mathrm{SHRINK}$	$\ll 1$ s	$\ll 1$ s
	LS1	≈ 40 s	≈ 1 : 20m
	Total	≈ 11 m	≈ 42 m
DIC	${\mathbf{x}}_{\mathrm{DIC}}$	≈ 1:24 h	≈ 7 m
	LS1	≈ 40 s	≈ 1:20 m
	Total	≈ 1:28 h	≈ 28 m
Proposed	${\mathbf{x}_{\mathrm{CNN}}}$	≈ 4 m	≈ 5 s
	LS1	≈ 40 s	≈ 1:20 m
	Total	≈ 8 m	≈ 21 m

Obviously, in terms of achieved image quality, the advantage of the DIC- and the CNN-based Tikhonov regularization are given by obtaining stronger priors which allow to use a smaller number of iterations to regularize the solution. The advantage of our proposed approach compared to the dictionary learning-based is the highly reduced time to compute the prior which is used for regularization. The reason lies in the fact that the DIC-based method requires to solve problem (23) to obtain the prior ${\mathbf{x}}_{\mathrm{DIC}}$, while in our method a CNN is used to obtain the prior ${\mathbf{x}_{\mathrm{CNN}}}$. Since problem (23) is separable, OMP is applied for each image/volume patch which is prohibitive as the number of overlapping patches in a 3D volume is in the order of $\mathcal{O}(N_x\cdot N_y \cdot N_z)$ or $\mathcal{O}(N_x\cdot N_y \cdot N_t)$, respectively. Obtaining ${\mathbf{x}_{\mathrm{CNN}}}$, on the other hand, does not involve the solution of any minimization problem but only requires the application of the network $u_{\theta}$ to the different patches. As this corresponds to matrix-vector multiplications with sparse matrices, its computational cost is lower and the calculations are further highly accelerated by performing the computations on a GPU.