Deep Convolutional Neural Network for Inverse Problems in Imaging

I Introduction

Over the past decades, iterative reconstruction methods have become the dominant approach to solving inverse problems in imaging including denoising , deconvolution , and interpolation . Thanks to robust regularizers such as total variation and sparsity , practical algorithms have appeared with excellent image quality and reasonable computational complexity. These advances have been particularly influential in the field of biomedical imaging, e.g., in magnetic resonance imaging (MRI) and X-ray computed tomography (CT) . These devices face an unfavorable trade-off between noise and acquisition time. Short acquisitions lead to severe degradations of image quality, while long acquisitions may cause motion artifacts, patient discomfort, or even patient harm in the case of radiation-based modalities. Iterative reconstruction with regularization provides a way to mitigate these problems in software, i.e. without developing new scanners. With the appearance of compressed sensing , our theoretical understanding of these approaches evolved further, and, in some applications, remarkable outcomes appeared with stable reproducibility .

A more recent trend is deep learning , which has arisen as a promising framework providing state-of-the-art performance for image classification and segmentation . Moreover, regression-type neural networks demonstrated impressive results on inverse problems with exact models such as signal denoising, deconvolution, and interpolation . Central to this resurgence of neural networks has been the convolutional neural network (CNN) architecture. Whereas the classic multilayer perceptron consists of layers that can perform arbitrary matrix multiplications on their input, the layers of a CNN are restricted to perform convolutions, greatly reducing the number of parameters which must be learned.

Researchers have begun to investigate the link between conventional approaches and deep learning networks . Gregor and LeCun explored the similarity between the ISTA algorithm and a shared layerwise neural network and demonstrated that several layer-wise neural networks act as a fast approximated sparse coder. In , a nonlinear diffusion reaction process based on the Perona-Malik process was proposed using deep convolutional learning; convolutional filters from diffusion terms were trained instead of using well-chosen filters like kernels for diffusion gradients, while the reaction terms were matched to the gradients of a data fidelity term. In , the authors focused on the relationship between $l_{0}$ penalized-least-squares methods and deep neural networks. In the context of a clustered dictionary model, they found that the non-shared layer-wise independent weights and activations of a deep neural network provide more performance gain than the layer-wise fixed parameters of an unfolded $l_{0}$ iterative hard thresholding method. The quantitative analysis relied on the restricted isometry property (RIP) condition from compressed sensing . Others have investigated learning optimal shrinkage operators for deep-layered neural networks .

Despite these works, practical and theoretical questions remain regarding the link between iterative reconstruction and CNNs. For example, in which problems can CNNs outperform traditional iterative reconstructions, and why? Where does this performance come from, and can the same gains be realized by learning aspects of the iterative process (e.g. the shrinkage)? Although began to address this connection, they only assumed that the filters learned in the Perona-Malik scheme are modified gradient kernels, with performance gains coming from the increased size of the filters.

In this paper, we explore the relationship between CNNs and iterative optimization methods for one specific class of inverse problems: those where the normal operator associated with the forward model ($H^{*}H$, where $H$ is the forward operator and $H^{*}$ is the adjoint operator) is a convolution. The class trivially includes denoising and deconvolution, but also includes MRI , X-ray CT , and diffraction tomography (DT). Based on this connection, we propose a method for solving these inverse problems by combining a fast, approximate solver with a CNN. We demonstrate the approach on low-view CT reconstruction, using filtered back projection (FBP) and a CNN that makes use of residual learning and multilevel learning . We use high-view FBP reconstructions for training, meaning that training is possible from real data (without oracle knowledge). We compare to a state-of-the art regularized iterative reconstruction and show promising results on both synthetic and real CT data. Especially, reconstructed images from the proposed network represented well complex textures which are important to diagnosis.

II Inverse Problems with Shift-Invariant Normal Operators

We begin our discussion by describing the class of inverse problems for which the normal operator is a convolution. As we will show, solving problems of this form iteratively requires repeated convolutions and point-wise nonlinearities, which suggests that CNNs may offer an alternative solution. The class is broad, encompassing at least denoising, deconvolution, and reconstruction of MRI, CT, and diffraction tomography images. The underlying convolutional structure is known for MRI and CT and has been exploited in the past for the design of fast algorithms (e.g. ). Here, we aim to give a general and concise definition to motivate our method. We go on to discuss the standard direct and iterative approaches to solving these problems.

A convolution, $H_{h}:L_{2}(\Omega)\to L_{2}(\Omega)$, is a linear operator such that $H_{h}f=\mathcal{F}^{*}M_{\hat{h}}\mathcal{F}f$, where $\mathcal{F}$ is the Fourier transform, $\hat{h}$ is the Fourier transform of $h$, and $M_{\hat{h}}$ is a multiplication.

A reversible change of variables, $\Phi_{\varphi}:L_{2}(\Omega_{1})\to L_{2}(\Omega_{2})$, is a linear operator such that $\Phi_{\varphi}f=f(\varphi(\cdot))$ for some $\phi:\Omega_{2}\to\Omega_{1}$ and such that its inverse, $\Phi_{\varphi}^{-1}=\Phi_{\varphi^{-1}}$ exists.

If $H_{h}$ is a convolution, then $H_{h}^{*}H_{h}$ is as well (because $\mathcal{F}^{*}M_{\hat{h}}^{*}\mathcal{F}\mathcal{F}^{*}M_{\hat{h}}\mathcal{F}=\mathcal{F}^{*}M_{|\hat{h}|^{2}}\mathcal{F}$), but this is true for a wider set of operators. Theorem 1 describes this set.

If there exists an isometry, $T$, a multiplication, $M_{m}$, and a change of variables, $\Phi_{\varphi}$, such that $H=TM_{m}\Phi_{\varphi}^{-1}\mathcal{F}$, then $H^{*}H$ is a convolution with $\hat{h}=|\det J_{\varphi}|M_{\Phi_{\varphi}|m|^{2}}$, where $J_{\varphi}$ is the Jacobian matrix of $\varphi$ and $M_{\Phi_{\varphi}|m|^{2}}$ is a suitable multiplication.

Given an operator, $H$, that satisfies the conditions of Theorem 1,

where (a) follows from the definitions of isometry and multiplication and (b) follows from the definition of a reversible change of variables. Thus, $H^{*}H$ is a convolution by Definition 2. ∎

where $\Phi_{\varphi}$ changes from Cartesian to polar coordinates (i.e. $\varphi^{-1}(\theta,r)=(r\cos\theta,r\sin\theta)$) and $T$ is the inverse Fourier transform with respect to $r$ (which is an isometry due to Parseval’s theorem). This maps a function, $f$, of space, $\bm{x}$, to its Fourier transform, $\hat{f}$, which is a function of frequency, $\bm{\omega}$. Then, it performs a change of variables, giving $\hat{f}_{\text{polar}}$, which is a function of a polar frequency variables, $(\theta,r)$. Finally, $T$ inverts the Fourier transform along $r$, resulting in a sinogram that is a function of $\theta$ and a polar space variable, $y$. Theorem 1 states that $R^{*}R$ is a convolution with $\hat{h}(\bm{\omega})=|\det J_{\varphi}(\bm{\omega})|=1/\|\bm{\omega}\|$, where, again, $\bm{\omega}$ is the frequency variable associated with the 2D Fourier transform, $\mathcal{F}$.

II-B Direct Inversion

Given a normal-convolutional operator, $H$, the inverse (or reconstruction) problem is to recover an image $f$ from its measurements $g=Hf$. The theory presented above suggests two methods of direct solutions to this problem. The first is to apply the inverse of the filter corresponding to $H^{*}H$ to the back projected measurements,

where $W_{h}$ is a convolution operator with $\hat{h}(\bm{\omega})=1/(\lvert\det J_{\varphi}\rvert\Phi_{\varphi}\lvert m(\bm{\omega})\rvert^{2})$. This is exactly equivalent to performing a deconvolution in the reconstruction space. The second is to invert the action of $H$ in the measurement domain before back projecting,

where $M_{h}$ is a multiplication operator with $h(\bm{\omega})=1/(|\det J_{\varphi}||m(\bm{\omega})|^{2})$. If $T$ is a Fourier transform, then this inverse is a filtering operation followed by a back projection; if $T$ is not, the operation remains filtering-like in the sense that it is diagonalizable in the transform domain associated with $T$. Note also that if $T$ is not a Fourier transform, then the variable $\bm{\omega}$ no longer refers to frequency. Given the their filter-like form, we refer to these direct inverses as filtered back projection (FBP) , a term borrowed from X-ray CT reconstruction.

Returning to the example of the continuous 2D X-ray transform, the first method would be to back project the measurements and then apply the filter with a 2D Fourier transform given by $\|\bm{\omega}\|$. The second approach would be to apply the filter with 1D Fourier transform given by $\omega$ to each angular measurement and then back project the result. In the continuous case, the methods are equivalent, but, in practice, the measurements are discrete and applying these involves some approximation. Then, which form is used affects the accuracy of the reconstruction (along with the runtime). This type of error can be mitigated by formulating the FBP to explicitly include the effects of sampling and interpolation (e.g., as in ). The larger problem is that the filter greatly amplifies noise, thus in practice some amount of smoothing is also applied.

II-C Iterative Inversion

This formulation does not admit a closed form solution, and, therefore, is typically solved iteratively. For example, the popular ISTA algorithm solves Eq. (5) with the iterate

where $\mathcal{S}_{\theta}$ is the soft-thresholding operator by value $\theta$ and $L\leq\operatorname*{eig}({\mathbf{W}}^{*}{\mathbf{H}}^{*}{\mathbf{H}}{\mathbf{W}})$ is the Lipschitz constant of a normal operator. When the forward model is normal-convolutional and when ${\mathbf{W}}$ is a convolution, the algorithm consists of iteratively filtering by ${\mathbf{I}}-(1/L){\mathbf{W}}^{*}{\mathbf{H}}^{*}{\mathbf{H}}{\mathbf{W}}$, adding a bias, $(1/L){\mathbf{W}}^{*}{\mathbf{H}}^{*}{\mathbf{y}}$, and applying a point-wise nonlinearity, $\mathcal{S}_{\theta}$. This is illustrated as a block diagram with unfolded iterates in Fig. 1 (b). Many other iterate methods for solving Eq. (5), including ADMM , FISTA, and SALSA, also rely on these basic building blocks.

III Proposed Method: FBPConvNet

The success of iterative methods consisting of filtering plus pointwise nonlinearities on normal-convolutional inverse problems suggests that CNNs may be a good fit for these problems as well. Based on this insight, we propose a new approach to these problems, which we call the FBPConvNet. The basic structure of the FBPConvNet algorithm is to apply the discretized FBP from Section II-B to the measurements and then use this as the input of a CNN which is trained to regress the FBP result to a suitable ground truth image. This approach applies in principle to all normal-convolutional inverse problems, but we have focused in this work on CT reconstruction. We now describe the method in detail.

While it would be possible to train a CNN to regress directly from the measurement domain to the reconstruction domain, performing the FBP first greatly simplifies the learning. The FBP encapsulates our knowledge about the physics of the inverse problem and also provides a warm start to the CNN. For example, in the case of CT reconstruction, if the sinogram is used as input, the CNN must encode a change between polar and Cartesian coordinates, which is completely avoided when the FBP is used as input. We stress again that, while the FBP is specific to CT, Section II-C shows that efficient, direct inversions are always available for normal-convolutional inverse problems.

III-B Deep Convolutional Neural Network Design

While we were inspired by the general form of the proximal update, (6), to apply a CNN to inverse problems of this form, our goal here is not to imitate iterative methods (e.g. by building a network that corresponds to an unrolled version of some iterative method), but rather to explore a state-of-the-art CNN architecture. To this end, we base our CNN on the U-net , which was originally designed for segmentation. There are several properties of this architecture that recommend it for our purposes.

Multilevel decomposition. The U-net employs a dyadic scale decomposition based on max pooling, so that the effective filter size in the middle layers is larger than that of the early and late layers. This is critical for our application because the filters corresponding to $H^{*}H$ (and its inverse) may have non-compact support, e.g. in CT. Thus, a CNN with a small, fixed filter size may not be able to effectively invert $H^{*}H$. This decomposition also has a nice analog to the use of multiresolution wavelets in iterative approaches.

Multichannel filtering. U-net employs multichannel filters, such that there are multiple feature maps at each layer. This is the standard approach in CNNs to increase the expressive power of the network . The multiple channels also have an analog in iterative methods: In the ISTA formulation (6), we can think of the wavelet coefficient vector ${\mathbf{a}}$ as being partitioned into different channels, with each channel corresponding to one wavelet subband . Or, in ADMM, the split variables can be viewed as channels. The CNN architecture greatly generalizes this by allowing filters to make arbitrary combinations of filters.

Residual learning. As a refinement of the original U-net, we add a skip connection between input and output, which means that the network actually learns the difference between input and output. This approach mitigates the vanishing gradient problem during training. This yields a noticeable increase in performance compared to the same network without the skip connection.

Implementation details. We made two additional modification to U-net. First, we use zero-padding so that the image size does not decrease after each convolution. Second, we replaced the last layer with a convolutional layer which reduces the 64 channels to a single output image. This is necessary because the original U-net architecture results in two channgels: foreground and background.

IV Experiments and results

We now describe our experimental setup and results. Though the FBPConvNet algorithm is general, we focus here on sparse-view X-ray CT reconstruction. We compare FBPConvNet to FBP alone and a state-of-the-art iterative reconstruction method . This method (which we will refer to as the TV method for brevity) solves a version of Eq. (5) with the popular TV regularization via ADMM. It exploits the convolutional structure of $H^{*}H$ by using FFT-based filtering in its iterates.

Our experiments proceed as follows: We begin with a full view sinogram (either synthetically generated or from real data). We compute its FBP (standard high quality reconstruction) and take this as the ground truth. We then compare the results of applying the TV method to the subsampled sinogram with the results of applying the FBPConvNet to the same. This type of sparse-view reconstruction is of particular interest for human imaging because, e.g., a twenty times reduction in the number of views corresponds to a twenty times reduction in the radiation dose received by the patient.

We used three datasets for evaluations of the proposed method. The first two are synthetic in that the sinograms are computed using a digital forward model, while the last comes from real experiments.

The ellipsoid dataset is a synthetic dataset that comprises 500 images of ellipses of random intensity, size, and location. Sinograms for this data are 729 pixels by 1,000 views and are created using the analytical expression for the X-ray transform of an ellipse. The Matlab function iradon is used for FBPs.

The biomedical dataset is a synthetic dataset that comprises 500 real in-vivo CT images from the Low-dose Grand challenge competition from database made by the Mayo clinic. Sinograms for this data are 729 pixels by 1,000 views and are created using the Matlab function radon. iradon is again used for FBPs.

The experimental dataset is a real CT dataset that comprises 377 sinograms collected from an experiment at the TOMCAT beam line of the Swiss Light Source at the Paul Scherrer Institute in Villigen, Switzerland. Each sinogram is 1493 pixels by 721 views and comes from one z-slice of a single rat brain. FBPs were computed using our own custom routine which closely matches the behavior of iradon while accommodating different sampling steps in the sinogram an reconstruction domains.

To make sparse-view FBP images in synthetic datasets, we uniformly subsampled the sinogram by factors of 7 and 20 corresponding to 143 and 50 views, respectively. For the real data, we subsampled by factors of 5 and 14 corresponding to 145 and 52 views.

IV-B Training Procedure

FBPConvNet. In case of synthetic data, the total number of training images is 475. The number of test images is 25. In the case of the biomedical dataset, the test data is chosen from a different subject than the training set. For the real data, the total number of training images is 327. The number of test images is 25. The test data are obtained from the last z-slices with the gap of 25 slices left between testing and training data. All images are scaled between 0 and 550.

The CNN part of the FBPConvNet is trained using pairs of low-view FBP images and full-view FBP images as input and output, respectively. Note that this training strategy means that the method is applicable to real CT reconstructions where we do not have access to an oracle reconstruction.

We use the MatConvNet toolbox (ver. 20) to implement the FBPConvNet training and evaluation, with a slight modification: We clip the computed gradients to a fixed range to prevent the divergence of the cost function . We use a Titan Black GPU graphic processor (NVIDIA Corporation) for training and evaluation. Total training time is about 15 hours for 101 iterations.

The hyper parameters for training are as follows: learning rate decreasing logarithmically from 0.01 to 0.001; batchsize equals 1; momentum equals 0.99; and the clipping value for gradient equals $10^{-2}$. We use flip-flop data augmentation in both horizontal and vertical directions during the training phase to reduce overfitting .

State-of-the-art TV reconstruction. For completeness, we comment on how the iterative method used the training and testing data. Though it may be a fairer comparison to require the TV method to select its parameters from the training data (as the FBPConvNet does), we instead simply select the parameters that optimize performance on the training set (via a golden-section search). We do this with the understanding that the parameter is usually tuned by hand in practice and because the correct way to learn these parameters from data remains an open question.

V Experimental Results

We use SNR as a quantitative metric. If ${\mathbf{x}}$ is the oracle and $\hat{{\mathbf{x}}}$ is the reconstructed image, SNR is given by

where a higher SNR value corresponds to a better reconstruction.

Figures 3 and 4 and Table I show the results for the ellisoidal dataset. In the seven times downsampling case, Figure 3, the full-view FBP (ground truth) shows nearly artifact-free ellipsoids, while the sparse-view FBP shows significant line artifacts (most visible in the background). Both the TV and FBPConvNet methods significantly reduce these artifacts, giving visually indistinguishable results. When the downsampling is increased to twenty times, Figure 4, the line artifacts in the sparse-view FBP reconstruction are even more pronounced. Both the TV and FBPConvNet reduce these artifacts, though the FBPConvNet retains some of the artifacts. The average SNR on the testing set for the TV method is higher than that of the the FBPConvNet. This is a reasonable results given that the phantom is piecewise constant and thus the TV regularization should be optimal .

V-B Biomedical Dataset

Figures 5 and 6 and Table II show the results for the biomedical dataset. In Figure 5, again, the sparse-view FBP contains line artifacts. Both TV and the proposed method remove streaking artifacts satisfactorily; however, the TV reconstruction shows the cartoon-like artifacts that are typical of TV reconstructions. This trend is also observed in severe case (x20) in Fig. 6. Quantitatively, the proposed method outperforms the TV method.

V-C Experimental Dataset

Figures 7 and 8 and Table III show the results for the experimental dataset. The SNRs of all methods are significantly lower here because of the relatively low contrast of the sinogram. In Fig. 7, we observe the same trend as for the biomedical dataset, where the TV method oversmooths and the FBPConvNet better preserves fine structures. These trends also appears in twenty times downsampling case (x20) in Fig. 8. The FBPConvNet had a higher SNR than the TV method in both settings.

VI Discussion

The experiments provide strong evidence for the feasibility of the FBPConvNet for sparse-view CT reconstruction. The conventional iterative algorithm with TV regularization outperformed the FBPConvNet in the ellipsoidal dataset, while the reverse was true for the biomedical and experimental datasets. In these more-realistic datasets, the SNR improvement of the FBPConvNet came from its ability to preserve fine details in the images. This points to one advantage of the proposed method over iterative methods: the iterative methods must explicitly impose regularization, while the FBPConvNet effectively learns a regularizer from the data.

The computation time for the FBPConvNet was about 200 ms for the FBP and 200$\sim$300 ms in GPU for the CNN for a $512\times 512$ image. This is much faster than the iterative reconstruction, which, in our case, requires around 7 minutes even after the regularization parameters have been selected.

A major limitation of the proposed method is lack of transfer between datasets. For instance, when we put FBP images from a twenty-times subsampled sinogram into the network trained on the seven-times subsampled sinogram, the results retain many artifacts. Handling datasets of different dimensions or subsampling factors requires retraining the network. Future work could address strategies for heterogeneous datasets.

Our theory suggests that the methodology proposed here is applicable to all problems where the normal operator is shift-invariant; however, we have focused here on X-ray CT reconstruction. We expect that adapting the method to, e.g., MRI reconstruction will be non-trivial experimentally, because it will require large sets of training data (either from a high-quality forward model or real data) and a high-quality iterative reconstruction algorithm for comparison. Furthermore, because MRI and DT involve complex values (both in the measurement and reconstruction domains), we need a CNN architecture that correctly handles complex values. Therefore, we leave experimentation on other modalities to future work.

VII Conclusion

In this paper, we proposed a deep convolutional network for inverse problems with a focus on biomedical imaging. The proposed method, which we call the FBPConvNet combines FBP with a multiresolution CNN. The structure of the CNN is based on U-net, with the addition of residual learning.

This approach was motivated by the convolutional structure of several biomedical inverse problems, including CT, MRI, and DT. Specifically, we showed conditions on a linear operator that ensure that its normal operator is a convolution. This results suggests that CNNs are well-suited to this subclass of inverse problems.

The proposed method demonstrated compelling results on synthetic and real data. It compared favorably to state-of-the-art iterative reconstruction on the two more realistic datasets. Furthermore, after training, the computation time of the proposed network per one image is under a second.

Acknowledgment

The authors would like to thank Dr. Cynthia McCollough, the Mayo Clinic, the American Association of Physicists in Medicine, and grants EB017095 and EB017185 from the National Institute of Biomedical Imaging and Bioengineering for giving opportunities to use real-invivo CT DICOM images (Fig. 5-6). The authors also thank thank Dr. Marco Stampanoni, Swiss Light Source, Paul Scherrer Institute, Villigen, Switzerland, for providing real CT sinograms (Fig. 7-8).

References