Unpaired image denoising using a generative adversarial network in X-ray CT

TL;DR

This paper introduces a novel unpaired deep learning approach using a fidelity-embedded GAN for denoising low-dose X-ray CT images, effectively preserving features without requiring paired training data.

Contribution

It develops a GAN-based denoising method that learns from unpaired low-dose and standard-dose CT images, enabling effective noise reduction while maintaining image details.

Findings

Method performs comparably to paired-data approaches.

Preserves fine anomalous features during denoising.

Validated with clinical experiments on low-dose CT images.

Abstract

This paper proposes a deep learning-based denoising method for noisy low-dose computerized tomography (CT) images in the absence of paired training data. The proposed method uses a fidelity-embedded generative adversarial network (GAN) to learn a denoising function from unpaired training data of low-dose CT (LDCT) and standard-dose CT (SDCT) images, where the denoising function is the optimal generator in the GAN framework. This paper analyzes the f-GAN objective to derive a suitable generator that is optimized by minimizing a weighted sum of two losses: the Kullback-Leibler divergence between an SDCT data distribution and a generated distribution, and the $\ell_2$ loss between the LDCT image and the corresponding generated images (or denoised image). The computed generator reflects the prior belief about SDCT data distribution through training. We observed that the proposed method allows the preservation of fine anomalous features while eliminating noise. The experimental results show that the proposed deep-learning method with unpaired datasets performs comparably to a method using paired datasets. A clinical experiment was also performed to show the validity of the proposed method for noise arising in the low-dose X-ray CT.

Tables

Table 1. Table 1: PERFORMANCE COMPARISON BETWEEN THE RESULTS FROM EXECUTING THE PROPOSED METHOD WITH GAUSSIAN and QUANTUM NOISE, USING PAIRED AND UNPAIRED DATASETS.

Measures	Guassian noise			Quantum noise
	LDCT	paired	unpaired	LDCT	paired	unpaired
PSNR(dB)	24.3305	35.2691	34.8404	23.9892	33.5025	33.4212
SSIM	0.2600	0.9154	0.9147	0.3176	0.8831	0.8688
MSE	1.3284e+03	109.4953	120.2068	1.4368e+03	161.0430	164.7475

Table 2. Table 2: Performance comparison of the proposed method for the test datasets acquired at different noise levels

	$I_0=5\ast {10}^4$		$I_0=6.5\ast {10}^4$		$I_0=8\ast {10}^4$		$I_0=9.5\ast {10}^4$
measure	LDCT	unpair	LDCT	unpair	LDCT	unpair	LDCT	unpair
PSNR	22.1787	27.9628	23.1996	31.3124	23.9892	33.4212	24.6277	33.8760
SSIM	0.2425	0.8639	0.2826	0.8756	0.3176	0.8688	0.3486	0.8884
MSE	2.1801e+03	578.6777	1.7233e+03	266.5740	1.4368e+03	164.7475	1.2404e+03	150.9356

Table 3. Table 3: Quantitative evaluation of the proposed method for clinical brain CT images.

	LDCT images	Corrected images	p-value
CT number (Mean$\pm$SD)
GM	32.45$\pm$1.22	31.78$\pm$1.06	<0.001
HM	26.36$\pm$1.52	26.34$\pm$1.08	0.81
SNR (Mean$\pm$SD)
GM	4.13$\pm$0.48	6.44$\pm$0.69	<0.001
HM	3.49 $\pm$0.41	5.53$\pm$0.59	<0.001
CNR (Mean$\pm$SD)
GM-HM	0.56 $\pm$0.19	0.78$\pm$0.23	<0.001

Equations

$$L_z(y):=\log p_x(y)-\lambda \Vert y-z{\Vert }_2^2$$

$$G_{\text{mbox}MAP}:=\arg \underset{G}{\min }{E}_{z\sim p_z}\left[-\log (p_x(G(z)))+\lambda ||G(z)-z|{|}^2\right].$$

$$G_{\text{mbox}MAP}$$

$$+\lambda E_{z\sim p_z}\left[||G(z)-z|{|}^2\right]$$

$$=\arg \underset{G}{\min }H(p_{G(z)},p_x)+\lambda E_{z\sim p_z}\left[||G(z)-z|{|}^2\right]$$

$$G_{\ast }:=\arg \underset{G}{\min }KL(p_{G(z)},p_x)+\lambda E_{z\sim p_z}\left[||G(z)-z|{|}^2\right]$$

$$G_{\ast }=\underset{G}{\mathrm{argmin}}\underset{D}{\max }J(D,G)$$

$$J(D,G):=E_{x\sim p_x}\left[D(x)\right]$$

$$+\lambda E_{z\sim p_z}\left[||G(z)-z|{|}^2\right]$$

$$D_{opt}=\arg \underset{D}{\max }{E}_{x\sim p_x}\left[D(x)\right]+E_{z\sim p_z}\left[\log (1-D(G(z)))\right]$$

$$D_{opt}$$

$$=\arg \underset{D}{\max }\int p_x(w)D(w)+p_{G(z)}(w)\log (1-D(w))dw$$

$$D_{opt}(w):=1-\frac{p_{G(z)}(w)}{p_x(w)},$$

$$\frac{d}{dt}\left(p_x(w)t+p_{G(z)}(w)\log (1-t)\right)=0.$$

$$E_{x\sim p_x}\left[D_{opt}(x)\right]+E_{z\sim p_z}\left[\log (1-D_{opt}(G(z)))\right]$$

$$+\lambda E_{z\sim p_z}\left[||G(z)-z|{|}^2\right]$$

$$=\int \left[p_x(w)D_{opt}(w)+p_{G(z)}(w)\log (1-D_{opt}(w))\right]dx$$

$$+\lambda E_{z\sim p_z}||G(z)-z|{|}^2$$

$$=\int \left[p_x(w)-p_{G(z)}(w)+p_{G(z)}(w)\log \left(\frac{p_{G(x)}(w)}{p_x(w)}\right)\right]dw$$

$$+\lambda E_{z\sim p_z}||G(z)-z|{|}^2$$

$$=KL(p_{G(z)},p_x)+\lambda E_{z\sim p_z}||G(z)-z|{|}^2,$$

$$G_{\ast \ast }=$$

$$\arg \underset{y}{\max }\log p_{x^{\ast }|z}(y|z)=\arg \underset{y}{\max }\left[\log p_x(y)+\log p_{z|x^{\ast }}(z|y)\right].$$

$$p_{z|x^{\ast }}(z|y)=p_{\epsilon |x^{\ast }}(z-y|y)=p_{\epsilon }(z-y)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{\Vert z-y{\Vert }_2^2}{2{\sigma }^2}}.$$

$$\arg \underset{y}{\max }\log p_{x^{\ast }|z}(y|z)$$

Full text

Computerized tomography, denoising, low-dose, generative adversarial network, unsupervised learning

I Introduction

In computed tomography (CT), reducing the X-ray radiation dose, while maintaining the diagnostic image quality is an important ongoing issue. This is because of growing concerns regarding the risk of radiation induced cancer Gonzalez2004 ; Miglioretti2013 . Low dose CT (LDCT) is commonly achieved by reducing the X-ray tube current. However, LDCT images obtained from commercial CT scanners in general suffer from the low signal-to-noise ratio (SNR) and the reduced diagnostic reliability. Therefore, numerous efforts have sought to denoise LDCT images, by finding a denoising mapping that converts an LDCT image to the corresponding standard-dose CT (SDCT) image.

Various iterative reconstruction (IR) methods have been proposed to reduce noise in LDCT images while preserving structure. The noise reduction modeling in these methods can employ loss functions in image space or in sinogram space. Some have incorporated prior knowledge into the denoising process by employing regularization strategies such as total variation (TV) Dong2014 ; Xu2011 , fractional-order TV Zhang2016 , and nonlocal TV Liu2016 . Markov random fields theory Li2004 ; Wang2006 or nonlocal means Chen2009 ; Ma2012 ; Zhang2017 have also been used as prior information. Statistical image reconstruction methods such as the maximum a posteriori (MAP) approach for data fitting in sinogram space have been used for efficient noise filtering of the sinogram in LDCT Sauer1993 ; Li2004 .

While these IR methods can significantly improve the quality of reconstructed CT images, they retain some limitations in clinical practice. To begin with, it is challenging to design a prior that conveys the characteristics of LDCT and SDCT images. For example, commonly used priors such as TV and its variants produce an over-smoothing effect that causes the loss of fine detail such as small anomalies. Next, IR methods impose high computational costs, as they require an iterative solver to find a reasonable approximate solution. Finally, sinogram-based methods use projection data for data fidelity, but it is generally difficult to access projection data from a commercial CT scanner.

Various supervised learning approaches have recently been suggested to reduce noise in LDCT images Chen2017a ; Chen2017 ; Kang2016 . Paired training data (i.e., LDCT and SDCT images) are not available in clinical practice, because it is in general difficult to obtain both types of images simultaneously from a given patient. Therefore, the learning approaches obtain paired training data by generating simulated LDCT image data from clinical SDCT images. Results with supervised datasets have shown that these approaches can reduce noise and artifacts in the simulated LDCT images. However, their practical performance depends heavily on the quality of the simulated LDCT image data. On the other hand, it is easy to collect a sufficient number of ”unpaired” SDCT and LDCT images from an in-hospital database, as the use of LDCT increases in routine clinical practice Chang2013 .

Motivated by the absence of paired data in the clinical environment, this study uses an unpaired learning approach to develop a noise reduction function that maps from LDCT images to its SDCT-like image counterparts. We exploit a method of MAP estimation for LDCT image denoising that imposes a compromise between noise reduction and a model-fitness with SDCT-like image priors obtained from the data of SDCT images. The problem of the MAP estimation can be converted into a minimization problem involving cross-entropy of data distribution and the distribution generated by the generative adversarial network (GAN) Goodfellow2014 . Given the general difficulty of directly handling cross-entropy, it is approximated to the Kullback-Leibler(KL) divergence, thus allowing an approximate MAP estimation using the GAN framework. This facilitates training of the network architecture with unpaired datasets.

The proposed GAN approach is capable of inferring the desired SDCT-like image priors from the data of the SDCT images. It is critical to choose effective training datasets to reflect the appropriate image priors, in order to preserve detailed features of the image while eliminating noise. This data dependency of SDCT image priors is validated by the experiment shown in Fig. 1. Under the assumption of Gaussian noise, the $\ell_{2}$ constraint is added to the proposed GAN model. Unlike conventional approaches that treat the entire image, the proposed method is carried out patch-by-patch manner, which can effectively train local noise features in the CT images. Numerical simulation and clinical results demonstrate that the proposed method has a great potential to reduce noise in LDCT images for unpaired CT images or even for a complicated noise model.

II Method

Throughout this paper, we denote by $z$ and $x$ LDCT and SDCT images, respectively, both of which are $n\times n$ pixel grayscale images. We assume that unpaired training samples $\{z_{k}\}_{k=1}^{N}$ and $\{x_{k}\}_{k=1}^{M}$ are drawn from unknown LDCT data distribution $p_{z}(z)$ and unknown SDCT data distribution $p_{x}(x)$, respectively. Here, $N$ and $M$ denote the number of LDCT and SDCT training samples, respectively. In terms of image denoising, $z\sim p_{z}$ corresponds to noisy image (source) and $x\sim p_{x}$ is regarded as noise-free image (target).

The object of our method for denoising is to learn a generator $G$ with input $z\sim p_{z}$ using GAN’s framework together with the unpaired training data, in such a way that the generator’s distribution $p_{G(z)}$ is an optimal approximation of the clean image distribution $p_{x}$ and the generators’s fidelity $\|z-G(z)\|_{2}$ is reasonably small for every $z\sim p_{z}$. Here, $\|z\|_{2}$ stands for the standard $\ell_{2}$-norm of $z$.

We adopt the additive Gaussian noise model that the noisy LDCT image $z$ is decomposed into a desired denoised image $x^{*}\sim p_{x}$ and additive Gaussian noise $\varepsilon$ of variance $\sigma$ at each pixel, i.e., $z=x^{*}+\varepsilon$ with $\varepsilon\sim\mathcal{N}(0,\sigma)$. One can estimate the clean image $x^{*}$ in terms of maximum a posteriori approach: Given $z\sim p_{z}$, find $\hat{x}$ that maximizes the conditional probability $p_{x^{*}|z}(\hat{x}|z)$ Sonderby2017 . This $\hat{x}$ is $\arg\max_{y}L_{z}(y)$, where

[TABLE]

and $\lambda=\frac{1}{2\sigma^{2}}$. For more detail, we refer to Appendix (V.1). However, this approach cannot be used directly due to the unknown distribution $p_{x}$. Instead of solving the denoising problem one-by-one for each $z\sim p_{z}$, we look for the denoising function $G:z\to\hat{x}$ by using GAN’s framework together with many samples $\{z_{k}\}_{k=1}^{N}$ and $\{x_{k}\}_{k=1}^{M}$.

An optimal denoising function can be derived by maximizing the expectation $E_{z\sim p_{z}}L_{z}(G(z))$ with respect to generator $G$:

[TABLE]

With the notion of the cross entropy $H(p,q)$, we have the following expression:

[TABLE]

However, it is still difficult to handle the cross-entropy $H(p_{G(z)},p_{x})$ only from the training samples even if we have a complete information of $G$. Instead, we approximate $G_{\mbox{\tiny MAP}}$ by

[TABLE]

where the cross-entropy $H(p_{G(z)},p_{x})$ in (II) is replaced by KL-divergence $KL(p_{G(z)},p_{x})=H(p_{G(z)},p_{x})-H(p_{G(z)})$. Interested readers may refer to the paperSonderby2017 for the effect of this replacement.

Now, we are ready to explain our fidelity-embedded GAN model. In GAN’s framework, the unpaired datasets are used to learn the generator $G$ that minimizes $KL(p_{G(z)},p_{x})$, while the discriminator $D$ tries to distinguish between generated $G(z)$ and $x\sim p_{x}$. The following theorem presents the proposed model which solves the optimal $G_{*}$ in (4) exactly when $G$ and $D$ reach their optimum.

Theorem II.1

The optimal $G_{*}$ in (4) is achieved by

[TABLE]

where

[TABLE]

Proof.   Given a fixed $G$, we first maximize $J(D,G)$ in (II.1) with respect to $D$. Since the third term in (II.1) is independent of $D$, it is enough to find

[TABLE]

Writing $x$ and $G(z)$ in (7) to the variable $w$, we have

[TABLE]

where the integral spans over the space of all images.

For a fixed $w$, the value of $D(w)=D_{opt}(w)$ that maximizes the integrand can be shown to be

[TABLE]

because $t=D_{opt}(w)$ satisfies

[TABLE]

Then, it follows from (II.1) and (9) that

[TABLE]

where $KL(p,q)$ is the KL-divergence given by $KL(p,q):=\int p(x)\log\left(\frac{p(x)}{q(x)}\right)dx$. Here, the last equality follows from the fact that $\int p_{x}(w)dw=\int p_{G(z)}(w)dw=1$. This completes the proof.

II.1 Data-driven dependency

The unpaired training samples $\mathcal{S}_{z}=\{z_{k}\}_{k=1}^{N}$ and $\mathcal{S}_{x}=\{x_{k}\}_{k=1}^{M}$ are used to compute the generator $G_{*}$ approximately:

[TABLE]

with $D=\arg\min_{D}\frac{1}{M}\sum_{x\in\mathcal{S}_{x}}D(x)+\frac{1}{N}$ $\sum_{z\in\mathcal{S}_{z}}\log(1-D(G(z)))$.

It should be noted that the computed generator $G_{**}$ depends heavily on the training samples. Fig. 1 shows the data-driven dependency, where we use three different training samples $\mathcal{S}^{j}=\{\mathcal{S}_{z},\mathcal{S}^{j}_{x}\},j=1,2,3$. The first training sample $\mathcal{S}^{1}$ consist of one disk and one rectangle with different sizes and positions (first column). The other samples ($\mathcal{S}^{2}$ and $\mathcal{S}^{3}$) are generated by adding a small anomaly to sample $\mathcal{S}^{1}$. In $\mathcal{S}^{2}$, the anomaly is added only in the rectangle ($\mathcal{S}^{2}$). On the other hand, $\mathcal{S}^{3}$ is generated by adding the anomaly only in the disk. In this case, the positions of the anomaly are randomly chosen. Fig. 1 shows the comparison of the performance of three generators $G_{**}^{1},~{}G_{**}^{2}$, and $G_{**}^{3}$, where $G_{**}^{j}$ is computed by (12) with $\mathcal{S}$ replaced by $\mathcal{S}^{j}$. The top rightmost image in Fig. 1 is a test image $z\sim p_{z}$, which contains an anomaly inside the disk. The $G_{**}^{1}(z)$ cannot preserve the anomaly inside the disk (the first row), whereas both $G_{**}^{2}(z)$ and $G_{**}^{3}(z)$ can preserve anomaly well. It is quite informative that $G_{**}^{2}(z)$ preserves the anomaly even though the training sample does not contain the anomaly inside the disk (the second row).

In this instance, the MAP estimation was completed in a patch-by patch manner, which is clearly explained in II.2.

II.2 Patch-based Network for Practical Applications

In practice, the proposed method is executed in a patch-by-patch manner, rather than working with the entire images (e.g., $512\times 512$ pixel CT images). Here, we take advantage of the fact Peyre2008 ; Peyre2009 that the patch manifolds of many images have a low dimensional structure. Moreover, the number of training data sets is significantly increased. These allow us to learn the generated distribution efficiently over the patches extracted from SDCT images. For simplicity, we will use the same notation of $x$ and $z$ to represent the image patches (there will be no confusion from the context).

Fig. 2 illustrates how the proposed method in (12) generates $G_{**}(z)$ from the unpaired training image patches. It optimizes $p_{G_{**}(z)}$ by forcing it to be close to $p_{x}$ while also minimizing the $\ell_{2}$ distance $\|G_{**}(z)-z\|_{2}$. Fig. 6 shows the performance of the $G_{**}$ with two different training samples; image patches of $\mathcal{S}^{1}$ and $\mathcal{S}^{2}$ are randomly selected, respectively, from 50 different CT slices (third column) and 600 different CT slices ($x$) (fourth column), respectively. As shown in Fig. 6, the performance of $G_{**}^{1}$ (using $\mathcal{S}^{1}$) is inferior to that of $G_{**}^{2}$ (using $\mathcal{S}^{2}$), owing to the lack of the diversity of $x$, which is called ‘mode collapse’ Che2017 ; Salimans2016 .

II.3 Network Architecture

For our generator, we adopt the deep convolutional framelet Ye2017 , which is a multi-scale convolutional neural network, consisting of a contracting path and an expansive path with skipped connections and a concatenation (concat) layer. Each step of the contracting and expansive paths consists of two repeated convolutions (conv) with a $3\times 3$ window. Each convolution is followed by a batch normalization (bnorm) and a leaky rectified linear unit (LReLU). Downsampling and upsampling of the features are performed by 2-D Haar wavelet decomposition (wave-dec) and recomposition (wave-rec), respectively. High pass filters after wavelet decomposition skip directly to the expansive path, while low pass filters (marked by ‘LF’ in Fig. 3) are concatenated with the features in the contracting path during the same step. At the end, an additional convolution layer is added to generate a grayscale output image. Note that each convolution in our network is performed with zero-padding to match the size of the input and output images. The architecture of the deep convolutional framelet is quite similar to that of the U-net Ronneberger2015 , a deep learning model widely used in image processing, except for the high pass filter pass. The deep convolutional framelet utilizes the wavelet decomposition and recomposition instead of pooling and unpooling operations for downsampling and upsampling, respectively. By doing so, more detailed information of image can be preserved during downsampling Ye2017 .

In an adversarial architecture, we adopt the PatchGAN classifier Isola2017 ; Li2016 ; Zhu2017 as a discriminator, which tries to classify whether each patch in an image is real or fake. PatchGAN enables us to learn the structural detail in images more easily than does the $1\times 1$ PixelGAN. The discriminator contains three convolution layers with a $4\times 4$ window and strides of two in each direction of the domain and each layer is followed by a batch normalization and a leaky ReLU with a slope of 0.2. As the final stage of the architecture, a $1\times 1$ convolution layer is added to generate 1-dimensional output data. Fig. 3 illustrates the architecture of the proposed method.

In our numerical simulation and clinical experiment, the objective function in (12) is minimized using an Adam optimizer Kingma2014 with a learning rate of 0.0002 and mini-batch size of 40, and 300 epochs are utilized for training. Training was implemented using Tensorflow Abadi2016 on a GPU (NVIDIA, Titan Xp. 12GB) system. It takes about a day to train the network. We empirically choose $\lambda=10$ in Eq. (12). The network weights were initialized following a Gaussian distribution with a mean of 0 and a standard deviation of 0.01.

II.4 Datasets

In our simulations and clinical experiments, all CT images of size $512\times 512$ are acquired by a 64-channel multi-detector CT scanner (Sensation 64; Siemens Healthcare, Forchheim, Germany).

II.4.1 Simulation study

We collected SDCT images of 62 patients acquired at a tube voltage of 120 kVp and a tube current-time product of 200 mAs. For each patient, we select 20 CT images that included the liver and portal vein. LDCT images are generated by adding noise to the SDCT images of the 62 patients. For the chosen SDCT images of 30 patients, corresponding LDCT images, and LDCT images of another 30 patients are used for the paired and unpaired training datasets, respectively. The LDCT images of the remaining two patients (not used for training) are used as a test dataset. See Fig. 5. From each image of size $512\times 512$, the patches of size $128\times 128$ are extracted with strides of 8 in each direction in the image domain. Among them, 40 patches are randomly selected, and they are used as training image patches.

In this study, we generate LDCT images with two different types of noise; Gaussian noise with a standard deviation of 25 (HU), and quantum noise based on the CT acquisition model. We generate a synthetic sinogram according to the Poisson quantum noise zeng2010 . In our simulations, the blank scan flux per detector was set to $8\times 10^{4}$, similar to Chen2017a . Electronic noise, modeled as Gaussian noise LaRiviere2005 , is also added to the synthetic sinogram. LDCT images are reconstructed by the filtered backprojection Bracewell1967 , where we only considered mono-chromatic and parallel X-ray beams.

II.4.2 Clinical study

For training our network, we collect 200 brain CT images of 10 patients, acquired at a tube voltage of 120 kVp and a fixed tube current-time product of 200 mAs. We further collect another 200 CT images at a tube voltage of 100 kVp and a reference tube current-time product of 190 mAs. We refer to the former and latter images as low-dose CT (LDCT) images, and standard-dose CT (SDCT) images, respectively. For LDCT images (obtained with a reference mAs of 190), the mAs varies from 156 mAs to 200 mAs along the length of the scan, due to the use of an automatic exposure compensation system. Fig. 9 shows unpaired brain CT images used for training. Using a low tube voltage and a tube current increases the image noise as shown in Fig. 9. Sixty LDCT images of three other patients who are not used in the training were collected for testing. As in the simulation, we generate unpaired patches of size $128\times 128$ from both LDCT and SDCT images.

III Results

The numerical simulation and the clinical experiment are performed in order to demonstrate the validity of the proposed approach for noise reduction in CT images.

Fig. 4 and Fig. 7 show the denoising results for liver CT images with simulated Gaussian and quantum noise, respectively. The second and third columns show the Results of the proposed method using paired and unpaired LDCT and SDCT datasets, respectively. The second and fourth rows show the zoomed regions of interest (ROI) of the rectangular regions of the images in the first and third rows, respectively. As shown in Fig. 4 and Fig. 7, both unpaired and paired learning methods reduce noise artifacts clearly, while preserving the morphological structures of the tissues. In Fig. 4, the fourth column shows the results of the model of (12) without the fidelity term $\|G_{**}(z)-z\|_{2}$. The method without the constraint generates plausible SDCT images, but tends to remove small details. Note the yellow arrow in the fourth column of Fig. 4.

Table 1 shows the comparison between the results from executing the proposed method with Gaussian and quantum noise, using paired and unpaired datasets. The average peak signal to noise ratio (PSNR), structure similarity (SSIM), and mean square error (MSE) between the reference and LDCT/corrected images are calculated for the test dataset (40 LDCT images of 2 patients). The performance of the proposed method with unpaired datasets shows comparable image quality (or structure similarity) to that with paired datasets in terms of PSNR, MSE and SSIM. Overall, compared with the quality of the LDCT images, the performance of the proposed method with Gaussian noise is slightly better than that with quantum noise. This is because fidelity in (12) is optimal for Gaussian noise.

We provide a performance evaluation of the proposed method for different noise levels (by using different blank scan flux $I_{0}$). For training, we use LDCT images with $I_{0}=8\times 10^{4}$. For testing, we use LDCT images with different noise levels (by using $I_{0}=5\times 10^{4},6.5\times 10^{4},8\times 10^{4}$, and $9.5\times 10^{4}$). We intentionally carried out training and testing with different noise levels to demonstrate that the proposed method effectively remove noises without causing a deterioration of the underlying morphological information. Fig. 8 and Table 2 show the results of the proposed method for different noise levels. In Table 2, the average PSNR, SSIM, and MSE between references and LDCT/corrected images are calculated for the test dataset (40 LDCT images of 2 patients). In terms of MSE, SSIM and PSNR, the proposed method reduces the noise in LDCT images and achieves the best performance at the lowest noise level ($I_{0}=9.5\times 10^{4}$). In the case of the highest noise level ($I_{0}=5\times 10^{4}$), overall HU values in the corrected images are less than the reference HU values, resulting in the highest MSE and lowest PSNR. On the other hand, the SSIM is similar at all noise levels. The robustness of the proposed method can be improved if LDCT images at different noise levels are used for training, as investigated in Chen2017a .

Fig. 10 shows the performance of the clinical experiment for brain CT images. The first row shows the brain LDCT images that are acquired at 178, 175, and 167 mAs. Severe noise artifacts occur due to the low radiation dose. The second row shows the denoising results of the proposed method based on unpaired LDCT and SDCT images that are obtained from different patients. Some examples of unpaired brain CT images used for training are illustrated in Fig. 9.

For quantitative analysis, we compute the SNR and CNR for test images (60 LDCT images of 3 patients) that were not used for training. The results are summarized in Table 3. In this study, the mean and SD of CT number are measured at manually selected ROIs on white matter (WM) and gray matter (GM) structures (Fig. 10). The size and shape of each ROI are kept constant. Because no reference image is available for each LDCT image, the SNR for each ROI is calculated as the mean attenuation (in Hounsfield units) divided by the SD Rapalino2012 . The CNR between the GM and WM regions is calculated as the difference in the mean attenuation of the two regions divided by the square root of the sum of their variances Rapalino2012 . A paired t-test with a significance level of $p<0.05$ is performed to assess the statistical significance between the LDCT images and the corrected images. The proposed method increases the SNR value of GM from $4.13\pm 0.48$ to $6.44\pm 0.69~{}(p<0.001)$, SNR value of HM from $3.49\pm 0.41$ to $5.53\pm 0.59~{}(p<0.001)$, and CNR value from $0.56\pm 0.19$ to $0.78\pm 0.23~{}(p<0.001)$. These results demonstrate that the proposed method reduces quantum noise in the corrected images. Even for general noise artifacts (not modeled by Gaussian noise), the proposed method significantly reduces the noise in terms of SNR and CNR.

IV Discussion and Conclusion

This study proposes an unpaired deep-learning method of image denoising for LDCT. The method approximates MAP estimation using a GAN framework, and can incorporate prior information on target SDCT images from a data distribution. Under the assumption of Gaussian noise, the $\ell_{2}$ constraint is added to the GAN framework. Numerical simulations show that the proposed deep learning method with unpaired datasets performs at least compatibly to supervised methods using paired datasets, in terms of PSNR and MSE (see Fig. 4). Clinical results also show that the proposed method enhances LDCT image quality even with unpaired SDCT and LDCT images (see Fig. 10). This type of approach was first used in amortized MAP inference Sonderby2017 to solve the super-resolution problem. Note that the proposed objective function (II.1) is different from that of Sonderby2017 as it is derived from the definition of $f$-GAN Nowozin2016 .

The proposed approach retains some issues requiring further research. One such issue is illustrated by the denoising effect with respect to the choice of samples from the data distribution in Figs. 1 and 6. This shows that, depending on the samples used for training, the proposed method can remove details, produce a plausible fake, or simply collapse when generating the denoised images. Therefore, it is essential for the image denoising technique to choose effective training datasets to reflect the appropriate image priors.

Note that the additive Gaussian noise model in (1) is not accurate enough. We adopted this model to simplify the mathematical analysis and to facilitate the derivation of the proposed method. Nevertheless, according to our experiments, the proposed method can learn the target SDCT image features from the data distribution, thus preserving the morphologically important information of the LDCT image and alleviating quantum noise, as well as Gaussian noise, as shown in Fig. 7 and 10. The proposed approach can be improved by applying a denoising algorithm that incorporates an accurate noise model of LDCT; for example, other constraints such as perceptual loss Yang2018 can be applied to the GAN framework in (12) to learn high-level features. In addition, an effective computational method is required for clinical applications. Given that the proposed method works a patch-by-patch, it is more time-consuming than conventional methods that work on an entire image at a time; it takes a few seconds to obtain a single $512\times 512$ denoised CT image. This can be reduced by using, for example, parallel computation.

Other than the issues mentioned above, future work will also focus on broadening the scope of the proposed method to provide a solution to the other CT reconstruction problems. Such problems include a reduction of the metal artifacts and scattering. This would grant the proposed method the potential to resolve the fundamental difficulty in applying deep learning approaches in X-ray CT (i.e., collecting paired CT image datasets). In addition, it would be interesting to focus on the quantitative analysis of the relationship between the estimated image and the training datasets.

V Appendix

V.1 Derivation of MAP for additive Gaussian noise

We follow the assumption of Section (II). That is, the noisy image $z$ is decomposed into a desired denoised image $x^{*}\sim p_{x}$ and additive Gaussian noise $\varepsilon$ of variance $\sigma$ at each pixel. For the MAP approach, we want to estimate $x^{*}$ by maximizing the conditional probability $p_{x^{*}|z}(y|z)$ with respect to the estimation $y$. Bayes rule provides that

[TABLE]

Note that $p_{x}(y)=p_{x^{*}}(y)$. Since $\varepsilon$ is independent of the original image $x^{*}$, we have

[TABLE]

It follows from (13) and (14) that the denoised image $y$ is

[TABLE]

where $\lambda=\frac{1}{2\sigma^{2}}$.