# Low Dose CT Image Denoising Using a Generative Adversarial Network with
Wasserstein Distance and Perceptual Loss

Qingsong Yang, Pingkun Yan*, , Yanbo Zhang, , Hengyong Yu, , Yongyi Shi,
Xuanqin Mou, , Mannudeep K. Kalra, Yi Zhang , Ling Sun, and Ge Wang Asterisk
indicates corresponding author.This work was supported in part by the National
Natural Science Foundation of China under Grant 61671312, and in part by the
National Institute of Biomedical Imaging and Bioengineering/National
Institutes of Health under Grant R01 EB016977 and Grant U01 EB017140.The
authors would also like to thank NVIDIA Corporation for the donation of the
Titan Xp GPU used for this research.Q. Yang, P. Yan* and G. Wang are with
Department of Biomedical Engineering, Rensselaer Polytechnic Institute, Troy,
NY, 12180 (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>wangg6@rpi.edu)Y. Zhang and
H. Yu are with the Department of Electrical and Computer Engineering,
University of Massachusetts Lowell, Lowell, MA 01854
<EMAIL_ADDRESS>hengyong-yu@ieee.org)Y. Shi and X. Mou are with
the Institute of Image Processing and Pattern Recognition, Xian Jiaotong
University, Xian, Shaanxi 710049, China (email: xqmou@mail.xjtu.edu.cn)M. K.
Kalra is with Department of Radiology, Massachusetts General Hospital, Harvard
Medical School, Boston, MA, USA(email<EMAIL_ADDRESS>Y. Zhang is
with the College of Computer Science, Sichuan Universiyt, Chengdu 610065,
China (e-mail: yzhang@scu.edu.cn)L. Sun is with Huaxi MR Research Center
(HMRRC), Department of Radiology, West China Hospital, Sichuan University,
Chengdu 610041, China (e-mail<EMAIL_ADDRESS>
###### Abstract

The continuous development and extensive use of CT in medical practice has
raised a public concern over the associated radiation dose to the patient.
Reducing the radiation dose may lead to increased noise and artifacts, which
can adversely affect the radiologist’s judgement and confidence. Hence,
advanced image reconstruction from low-dose CT data is needed to improve the
diagnostic performance, which is a challenging problem due to its ill-posed
nature. Over the past years, various low-dose CT methods have produced
impressive results. However, most of the algorithms developed for this
application, including the recently popularized deep learning techniques, aim
for minimizing the mean-squared-error (MSE) between a denoised CT image and
the ground truth under generic penalties. Although the peak signal-to-noise
ratio (PSNR) is improved, MSE- or weighted-MSE-based methods can compromise
the visibility of important structural details after aggressive denoising.
This paper introduces a new CT image denoising method based on the generative
adversarial network (GAN) with Wasserstein distance and perceptual similarity.
The Wasserstein distance is a key concept of the optimal transport theory, and
promises to improve the performance of GAN. The perceptual loss suppresses
noise by comparing the perceptual features of a denoised output against those
of the ground truth in an established feature space, while the GAN focuses
more on migrating the data noise distribution from strong to weak
statistically. Therefore, our proposed method transfers our knowledge of
visual perception to the image denoising task and is capable of not only
reducing the image noise level but also trying to keep the critical
information at the same time. Promising results have been obtained in our
experiments with clinical CT images.

###### Index Terms:

Low dose CT, Image denoising, Deep learning, Perceptual loss, WGAN

##  I Introduction

X-ray computed tomography (CT) is one of the most important imaging modalities
in modern hospitals and clinics. However, there is a potential radiation risk
to the patient, since x-rays could cause genetic damage and induce cancer in a
probability related to the radiation dose [1, 2]. Lowering the radiation dose
increases the noise and artifacts in reconstructed images, which can
compromise diagnostic information. Hence, extensive efforts have been made to
design better image reconstruction or image processing methods for low-dose CT
(LDCT). These methods generally fall into three categories: (a) sinogram
filtration before reconstruction [3, 4, 5], (b) iterative reconstruction [6,
7], and (c) image post-processing after reconstruction [8, 9, 10].

Over the past decade, researchers were dedicated to developing new iterative
algorithms (IR) for LDCT image reconstruction. Generally, those algorithms
optimize an objective function that incorporates an accurate system model [11,
12], a statistical noise model [13, 14, 15] and prior information in the image
domain. Popular image priors include total variation (TV) and its variants
[16, 17, 18], as well as dictionary learning [19, 20]. These iterative
reconstruction algorithms greatly improved image quality but they may still
lose some details and suffer from remaining artifacts. Also, they require a
high computational cost, which is a bottleneck in practical applications.

On the other hand, sinogram pre-filtration and image post-processing are
computationally efficient compared to iterative reconstruction. Noise
characteristic was well modeled in the sinogram domain for sinogram-domain
filtration. However, sinogram data of commercial scanners are not readily
available to users, and these methods may suffer from resolution loss and edge
blurring. Sinogram data need to be carefully processed, otherwise artifacts
may be induced in the reconstructed images.

Differently from sinogram denoising, image post-processing directly operates
on an image. Many efforts were made• in the image domain to reduce LDCT noise
and suppress artifacts. For example, the non-local means (NLM) method was
adapted for CT image denoising [8]. Inspired by compressed sensing methods, an
adapted K-SVD method was proposed [9] to reduce artifacts in CT images. The
block-matching 3D (BM3D) algorithm was used for image restoration in several
CT imaging tasks [10, 21]. With such image post-processing, image quality
improvement was clear but over-smoothing and/or residual errors were often
observed in the processed images. These issues are difficult to address, given
the non-uniform distribution of CT image noise.

The recent explosive development of deep neural networks suggests new thinking
and huge potential for the medical imaging field [22, 23]. As an example, the
LDCT denoising problem can be solved using deep learning techniques.
Specifically, the convolutional neural network (CNN) for image super-
resolution [24] was recently adapted for low-dose CT image denoising [25],
with a significant performance gain. Then, more complex networks were proposed
to handle the LDCT denoising problem such as the RED-CNN in [26] and the
wavelet network in [27]. The wavelet network adopted the shortcut connections
introducted by the U-net [28] directly and the RED-CNN [27] replaced the
pooling/unpooling layers of U-net with convolution/deconvolution pairs.

Despite the impressive denoising results with these innovative network
structures, they fall into a category of an end-to-end network that typically
uses the mean squared error (MSE) between the network output and the ground
truth as the loss function. As revealed by the recent work [29, 30], this per-
pixel MSE is often associated with over-smoothed edges and loss of details. As
an algorithm tries to minimize per-pixel MSE, it overlooks subtle image
textures/signatures critical for human perception. It is reasonable to assume
that CT images distribute over some manifolds. From that point of view, the
MSE based approach tends to take the mean of high-resolution patches using the
Euclidean distance rather than the geodesic distance. Therefore, in addition
to the blurring effect, artifacts are also possible such as non-uniform
biases.

To tackle the above problems, here we propose to use a generative adversarial
network (WGAN) [31] with the Wasserstein distance as the discrepancy measure
between distributions and a perceptual loss that computes the difference
between images in an established feature space [29, 30].

The use of WGAN is to encourage that denoised CT images share the same
distribution as that of normal dose CT (NDCT) images. In the GAN framework, a
generative network $G$ and a discriminator network $D$ are coupled tightly and
trained simultaneously. While the $G$ network is trained to produce realistic
images $G(\bm{z})$ from a random vector $\bm{z}$, the $D$ network is trained
to discriminate between real and generated images [32, 33]. GANs have been
used in many applications such as single image super-resolution [29], art
creation [34, 35], and image transformation [36]. In the field of medical
imaging, Nie et al. [37] proposed to use GAN to estimate CT image from its
corresponding MR image. Wolterink et al. [38] are the first to apply GAN
network for cardiac CT image denoising. And Yu et al. [39] used GAN network to
handle the de-alising problem for fast CS-MRI. Promising results were achieved
in these works. We will discuss and compare the results of those two networks
in Section III since the proposed network is closely related with their works.

Despite its success in these areas, GANs still suffer from a remarkable
difficulty in training [33, 40]. In the original GAN [32], $D$ and $G$ are
trained by solving the following minimax problem

$\min_{G}\max_{D}L_{\mathrm{GAN}}(D,G)=\mathbb{E}_{\bm{x}\sim P_{r}}[\log
D(\bm{x})]\\\ +\mathbb{E}_{\bm{z}\sim P_{z}}[\log\left(1-D(G(\bm{z}))\right)]$
(1)

where $\mathbb{E}(\cdot)$ denotes the expectation operator; $P_{r}$ and
$P_{z}$ are the real data distribution and the noisy data distribution. The
generator $G$ transforms a noisy sample to mimic a real sample, which defines
a data distribution, denoted by $P_{g}$. When $D$ is trained to become an
optimal discriminator for a fixed $G$, the minimization search for $G$ is
equivalent to minimizing the Jensen-Shannon (JS) divergence of $P_{r}$ and
$P_{g}$, which will lead to vanished gradient on the generator $G$ [40] and
$G$ will stop updating as the training continues.

Consequently, Arjovsky et al. [31] proposed to use the Earth-Mover (EM)
distance or Wasserstein metric between the generated image samples and real
data for GAN, which is referred to as WGAN, because the EM distance is
continuous and differentiable almost everywhere under some mild assumptions
while neither KL nor JS divergence is. After that, an improved WGAN with
gradient penalty was proposed [41] to accelerate the convergence.

The rationale behind the perceptual loss is two-fold. First, when a person
compares two images, the perception is not performed pixel-by-pixel. Human
vision actually extracts and compares features from images [42]. Therefore,
instead of using pixel-wise MSE, we employ another pre-trained deep CNN (the
famous VGG [43]) for feature extraction and compare the denoised output
against the ground truth in terms of the extracted features. Second, from a
mathematical point of view, CT images are not uniformly distributed in a high-
dimensional Euclidean space. They reside more likely in a low-dimensional
manifold. With MSE, we are not measuring the intrinsic similarity between the
images, but just their superficial differences in the brute-force Euclidean
distance. By comparing images according their intrinsic structures, we should
project them onto a manifold and calculate the geodesic distance instead.
Therefore, the use of the perceptual loss for WGAN should facilitate producing
results with not only lower noise but also sharper details.

In particular, we treat the LDCT denoising problem as a transformation from
LDCT to NDCT images. WGAN provides a good distance estimation between the
denoised LDCT and NDCT image distributions. Meanwhile, the VGG-based
perceptual loss tends to keep the image content after denoising. The rest of
this paper is organized as follows. The proposed method is described in
Section II. The experiments and results are presented in Section III. Finally,
relevant issues are discussed and a conclusion is drawn in Section IV.

##  II Methods

###  II-A Noise Reduction Model

Let $\bm{z}\in\mathbb{R}^{N\times N}$ denote a LDCT image and
$\bm{x}\in\mathbb{R}^{N\times N}$ denote the corresponding NDCT image. The
goal of the denoising process is to seek a function $G$ that maps LDCT
$\bm{z}$ to NDCT $\bm{x}$:

$G:\bm{z}\to\bm{x}$ (2)

On the other hand, we can also take $\bm{z}$ as a sample from the LDCT image
distribution $P_{L}$ and $\bm{x}$ from the NDCT distribution or the real
distribution $P_{r}$. The denoising function $G$ maps samples from $P_{L}$
into a certain distribution $P_{g}$. By varying the function $G$, we aim to
change $P_{g}$ to make it close to $P_{r}$. In this way, we treat the
denoising operator as moving one data distribution to another.

Typically, noise in x-ray photon measurements can be simply modeled as the
combination of Poisson quantum noise and Gaussian electronic noise. On the
contrary, in the reconstructed images, the noise model is usually complicated
and non-uniformly distributed across the whole image. Thus there is no clear
clue that indicates how data distributions of NDCT and LDCT images are related
to each other, which makes it difficult to denoise LDCT images using
traditional methods. However, this uncertainty of noise model can be ignored
in deep learning denoising because a deep neural network itself can
efficiently learn high-level features and a representation of data
distribution from modest sized image patches through a neural network.

###  II-B WGAN

Compared to the original GAN network, WGAN uses the Wasserstein distance
instead of the JS divergence to compare data distributions. It solves the
following minimax problem to obtain both $D$ and $G$ [41]:

$\min_{G}\max_{D}L_{\mathrm{WGAN}}(D,G)=-\mathbb{E}_{\bm{x}}[D(\bm{x})]+\mathbb{E}_{\bm{z}}[D(G(\bm{z}))]\\\
+\lambda\mathbb{E}_{\hat{\bm{x}}}[(||\nabla_{\hat{\bm{x}}}D(\hat{\bm{x}})||_{2}-1)^{2}],$
(3)

where the first two terms perform a Wasserstein distance estimation; the last
term is the gradient penalty term for network regularization; $\hat{\bm{x}}$
is uniformly sampled along straight lines connecting pairs of generated and
real samples; and $\lambda$ is a constant weighting parameter. Compared to the
original GAN, WGAN removes the $\log$ function in the losses and also drops
the last sigmoid layer in the implementation of the discriminator $D$.
Specifically, the networks $D$ and $G$ are trained alternatively by fixing one
and updating the other.

###  II-C Perceptual Loss

While the WGAN network encourages that the generator transforms the data
distribution from high noise to a low noise version, another part of the loss
function is added for the network to keep image details or information
content. Typically, a mean squared error (MSE) loss function is used, which
tries to minimize the pixel-wise error between a denoised patch $G(\bm{z})$
and a NDCT image patch $\bm{x}$ as [25, 26]

$L_{\mathrm{MSE}}(G)=\mathbb{E}_{(\bm{x},\bm{z})}\left[\frac{1}{N^{2}}||G\left(\bm{z}\right)-\bm{x}||_{F}^{2}\right],$
(4)

where $||\cdot||_{F}$ denotes the Frobenius norm. However, the MSE loss can
potentially generate blurry images and cause the distortion or loss of
details. Thus, instead of using a MSE measure, we apply a perceptual loss
function defined in a feature space

$L_{\mathrm{Perceptual}}(G)=\mathbb{E}_{(\bm{x},\bm{z})}\left[\frac{1}{whd}||\phi(G(\bm{z}))-\phi(\bm{x})||_{F}^{2}\right],$
(5)

where $\phi$ is a feature extractor, and $w$, $h$, and $d$ stand for the
width, height and depth of the feature space, respectively. In our
implementation, we adopt the well-known pre-trained VGG-19 network [43] as the
feature extractor. Since the pre-trained VGG network takes color images as
input while CT images are in grayscale, we duplicated the CT images to make
RGB channels before they are fed into the VGG network. The VGG-19 network
contains 16 convolutional layers followed by 3 fully-connected layers. The
output of the 16th convolutional layer is the feature extracted by the VGG
network and used in the perceptual loss function,

$L_{\mathrm{VGG}}(G)=\mathbb{E}_{(\bm{x},\bm{z})}\left[\frac{1}{whd}||VGG(G(\bm{z}))-VGG(\bm{x})||_{F}^{2}\right]$
(6)

For convenience, we call the perceptual loss computed by VGG network VGG loss.

Combining Eqs. (3) and (6) together, we get the overall joint loss function
expressed as

$\min_{G}\max_{D}L_{\mathrm{WGAN}}(D,G)+\lambda_{1}L_{\mathrm{VGG}}(G)$ (7)

where $\lambda_{1}$ is a weighting parameter to control the trade-off between
the WGAN adversarial loss and the VGG perceptual loss.

###  II-D Network Structures

The overall view of the proposed network structure is shown in Fig. 1. For
convenience, we name this network WGAN-VGG. It consists three parts. The first
part is the generator $G$, which is a convolutional neural network (CNN) of 8
convolutional layers. Following the common practice in the deep learning
community [44], small $3\times 3$ kernels were used in each convolutional
layer. Due to the stacking structure, such a network can cover a large enough
receptive field efficiently. Each of the first 7 hidden layers of $G$ have 32
filters. The last layer generates only one feature map with a single $3\times
3$ filter, which is also the output of $G$. We use Rectified Linear Unit
(ReLU) as the activation function.

The second part of the network is the perceptual loss calculator, which is
realized by the pre-trained VGG network [43]. A denoised output image
$G(\bm{z})$ from the generator $G$ and the ground truth image $\bm{x}$ are fed
into the pre-trained VGG network for feature extraction. Then, the objective
loss is computed using the extracted features from a specified layer according
to Eq. (6). The reconstruction error is then back-propagated to update the
weights of $G$ only, while keeping the VGG parameters intact.

The third part of the network is the discriminator $D$. As shown in Fig. 2,
$D$ has 6 convolutional layers with the structure inspired by others’ work
[43, 29, 30]. The first two convolutional layers have 64 filters, then
followed by two convolutional layers of 128 filters, and the last two
convolutional layers have 256 filters. Following the same logic as in $G$, all
the convolutional layers in $D$ have a small $3\times 3$ kernel size. After
the six convolutional layers, there are two fully-connected layers, of which
the first has 1024 outputs and the other has a single output. Following the
practice in [31], there is no sigmoid cross entropy layer at the end of $D$.

The network is trained using image patches and applied on entire images. The
details are provided in Section III on experiments.

Figure 1: The overall structure of the proposed WGAN-VGG network. In Part 1,
$n$ stands for the number of convolutional kernels and $s$ for convolutional
stride. So, $n32s1$ means the convolutional layer has 32 kernels with stride
1. Figure 2: The structure of the discriminator network. $n$ and $s$ have the
same meaning as in Fig. 1

0: Set hyper-parameters, $\lambda=10,\alpha=1\times
10^{-5},\beta_{1}=0.5,\beta_{2}=0.9,\lambda_{1}=0.1,\lambda_{2}=0.1$,

0: Set the number of total epochs, $N_{epoch}=100$, the number of iteration
for discriminator training, $N_{D}=4$, the batch size $m=128$, and image patch
size of $80\times 80$.

0: Initial discriminator parameters $w_{0}$, initial generator parameters
$\theta_{0}$

0: Load VGG-19 network parameters

1: for $num\\_epoch=0,...,N_{epoch}$ do

2: for $t=1,...,N_{D}$ do

3: Sample a batch of NDCT image patches $\\{\bm{x}^{(i)}\\}_{i=1}^{m}$, latent
LDCT patches $\\{\bm{z}^{(i)}\\}_{i=1}^{m}$, and random numbers
$\\{\epsilon^{(i)}\\}_{i=1}^{m}\sim\mathrm{Uniform}[0,1]$

4: for $i=1,...,m$ do

5:
$\hat{\bm{x}}^{(i)}\leftarrow\epsilon^{(i)}\bm{x}^{(i)}+(1-\epsilon^{(i)})G(\bm{z}^{(i)})$

6: $L^{(i)}(D)\leftarrow D(G(\bm{z}^{(i)}))-D(\bm{x}^{(i)})+\lambda(||\nabla
D(\hat{\bm{x}}^{(i)})||_{2}-1)^{2}$

7: end for

8: end for

9: Update $D$:
$w\leftarrow\text{Adam}(\nabla_{w}\frac{1}{m}\sum_{i=1}^{m}L^{(i)}(D),w,\alpha,\beta_{1},\beta_{2})$

10: Sample a batch of LDCT patches $\\{\bm{z}^{(i)}\\}_{i=1}^{m}$ and
corresponding NDCT patches $\\{\bm{x}^{(i)}\\}_{i=1}^{m}$,

11: for $i=1,...,m$ do

12:
$L^{(i)}(G)\leftarrow\lambda_{1}L_{\mathrm{VGG}}(\bm{z}^{(i)},\bm{x}^{(i)})-D(G(\bm{z}^{(i)}))$

13: end for

14: Update $G$,
$\theta\leftarrow\text{Adam}(\nabla_{\theta}\frac{1}{m}\sum_{i=1}^{m}L^{(i)}(G),w,\alpha,\beta_{1},\beta_{2})$

15: end for

Figure 3: Optimization procedure of WGAN-VGG network.

###  II-E Other Networks

For comparison, we also trained four other networks.

  * •

CNN-MSE with only MSE loss

  * •

CNN-VGG with only VGG loss

  * •

WGAN-MSE with MSE loss in the WGAN framework

  * •

WGAN with no other additive losses

  * •

Original GAN

All the trained networks are summarized in Table. I.

TABLE I: Summary of all trained networks: their loss functions and trainable networks. Network | Loss
---|---
CNN-MSE | $\min_{G}L_{\mathrm{MSE}}(G)$
WGAN-MSE | $\min_{G}\max_{G}L_{\mathrm{WGAN}}(G,D)+\lambda_{2}L_{\mathrm{MSE}}(G)$
CNN-VGG | $\min_{G}L_{\mathrm{VGG}}(G)$
WGAN-VGG | $\min_{G}\max_{G}L_{\mathrm{WGAN}}(G,D)+\lambda_{1}L_{\mathrm{VGG}}(G)$
WGAN | $\min_{G}\max_{G}L_{\mathrm{WGAN}}(G,D)$
GAN | $\min_{G}\max_{G}L_{\mathrm{GAN}}(G,D)$

(a)

(b)

(c)

Figure 4: Plots of validation loss versus the number of epochs during the
training of the 5 networks. (a) MSE loss convergence, (b) VGG loss convergence
and (c) Wasserstein estimation convergence.

##  III Experiments

###  III-A Experimental Datasets

We used a real clinical dataset authorized for “the 2016 NIH-AAPM-Mayo Clinic
Low Dose CT Grand Challenge” by Mayo Clinic for the training and evaluation of
the proposed networks [45]. The dataset contains 10 anonymous patients’
normal-dose abdominal CT images and simulated quarter-dose CT images. In our
experiments, we randomly extracted 100,096 pairs of image patches from 4,000
CT images as our training inputs and labels. The patch size is $64\times 64$.
Also, we extracted 5,056 pairs of patches from another 2,000 images for
validation. When choosing the image patches, we excluded image patches that
were mostly air. For comparison, we implemented a state-of-the-art 3D
dictionary learning reconstruction technique as a representative IR algorithm
[19, 20]. The dictionary learning reconstruction was performed from the LDCT
projection data provided by Mayo Clinic.

###  III-B Network Training

In our experiments, all the networks were optimized using Adam algorithm [46].
The optimization procedure for WGAN-VGG network is shown in Fig. 3. The mini-
batch size was 128. The hyper-parameters for Adam were set as
$\alpha=1e-5,\beta_{1}=0.5,\beta_{2}=0.9$, and we chose $\lambda=10$ as
suggested in [41], $\lambda_{1}=0.1,\lambda_{2}=0.1$ according to our
experimental experience. The optimization processes for WGAN-MSE and WGAN are
similar except that line 12 was changed to the corresponding loss function,
and for CNN-MSE and CNN-VGG, lines 2-10 were removed and line 12 was changed
according to their loss functions.

The networks were implemented in Python with the Tensorflow library[47]. A
NVIDIA Titan XP GPU was used in this study.

###  III-C Network Convergence

To visualize the convergence of the networks, we calculated the MSE loss and
VGG loss over the 5,056 image patches for validation according to Eqs. (4) and
(6) after each epoch. Fig. 4 shows the averaged MSE and VGG losses
respectively versus the number of epochs for the five networks. Even though
these two loss functions were not used at the same time for a given network,
we still want to see how their values change during the training. In the two
figures, both the MSE and VGG losses decreased initially, which indicates that
the two metrics are positively correlated. However, the loss values of the
networks in terms of MSE are increasing in the following order, CNN-
MSE$<$WGAN-MSE$<$WGAN-VGG$<$CNN-VGG (Fig. 4a), yet the VGG loss are in the
opposite order (Fig. 4b). The MSE and VGG losses of GAN network are
oscillating in the converging process. WGAN-VGG and CNN-VGG have very close
VGG loss values, while their MSE losses are quite different. On the other
hand, WGAN perturbed the convergence as measured by MSE but smoothly converged
in terms of VGG loss. These observations suggest that the two metrics have
different focuses when being used by the networks. The difference between MSE
and VGG losses will be further revealed in the output images of the
generators.

In order to show the convergence of WGAN part, we plotted the estimated
Wasserstein values defined as
$|-\mathbb{E}[D(\bm{x})]+\mathbb{E}[D(G(\bm{z}))]|$ in Eq. (3). It can be
observed in Fig. 4(c) that increasing the number of epochs did reduce the
W-distance, although the decay rate becomes smaller. For the WGAN-VGG curve,
the introduction of VGG loss has helped to improve the perception/visibility
at a cost of a compromised loss measure. For the WGAN and WGAN-MSE curves, we
would like to note that what we computed is a surrogate for the W-distance
which has not been normalized by the total number of pixels, and if we had
done such a normalization the curves would have gone down closely to zero
after 100 epochs.

(a) Full Dose FBP

(b) Quarter Dose FBP

(c) DictRecon

(d) GAN

(e) CNN-MSE

(f) CNN-VGG

(g) WGAN

(h) WGAN-MSE

(i) WGAN-VGG

Figure 5: Transverse CT images of the abdomen demonstrate a low attenuation
liver lesion (in the red box) and a cystic lesion in the upper pole of the
left kidney (in the blue box). This display window is [-160, 240]HU.

(a) Full Dose FBP

(b) Quarter Dose FBP

(c) DictRecon

(d) GAN

(e) CNN-MSE

(f) CNN-VGG

(g) WGAN

(h) WGAN-MSE

(i) WGAN-VGG

Figure 6: Zoomed ROI of the red rectangle in Fig. 5. The low attenuation liver
lesion with in the dashed circle represents metastasis. The lesion is
difficult to assess on quarter dose FBP recon (b) due to high noise content.
This display window is [-160, 240]HU.

(a) Full Dose FBP

(b) Quarter Dose FBP

(c) DictRecon

(d) GAN

(e) CNN-MSE

(f) CNN-VGG

(g) WGAN

(h) WGAN-MSE

(i) WGAN-VGG

Figure 7: Transverse CT images of the abdomen demonstrate small low
attenuation liver lesions. The display window is [-160, 240]HU.

(a) Full Dose FBP

(b) Quarter Dose FBP

(c) DictRecon

(d) GAN

(e) CNN-MSE

(f) CNN-VGG

(g) WGAN

(h) WGAN-MSE

(i) WGAN-VGG

Figure 8: Zoomed ROI of the red rectangle in Fig. 7 demonstrates the two
attenuation liver lesions in the red and blue circles. The display window is
[-160, 240]HU.

###  III-D Denoising Results

To show the denoising effect of the selected networks, we took two
representative slices as shown in Figs. 5 and 7. And Figs. 6 and 8 are the
zoomed regions-of-interest (ROIs) marked by the red rectangles in Figs. 5 and
7. All the networks demonstrated certain denoising capabilities. However, CNN-
MSE blurred the images and introduced waxy artifacts as expected, which are
easily observed in the zoomed ROIs in Figs. 6e and 8e. WGAN-MSE was able to
improve the result of CNN-MSE by avoiding over-smooth but minor streak
artifacts can still be observed especially compared to CNN-VGG and WGAN-VGG.
Meanwhile, using WGAN or GAN alone generated stronger noise (Figs. 6g and 8g)
than the other networks enhanced a few white structures in the WGAN/GAN
generated images, which are originated from the low dose streak artifact in
LDCT images, while on the contrary the CNN-VGG and WGAN-VGG images are
visually more similar to the NDCT images. This is because the VGG loss used in
CNN-VGG and WGAN-VGG is computed in a feature space that is trained previously
on a very large natural image dataset [48]. By using VGG loss, we transferred
the knowledge of human perception that is embedded in VGG network to CT image
quality evaluation. The performance of using WGAN or GAN alone is not
acceptable because it only maps the data distribution from LDCT to NDCT but
does not guarantee the image content correspondence. As for the lesion
detection in these two slices, all the networks enhance the lesion visibility
compared to the original noisy low dose FBP images as noise is reduced by the
different approaches.

As for iterative reconstruction technique, the reconstruction results depend
greatly on the choices of the regularization parameters. The implemented
dictionary learning reconstruction (DictRecon) result gave the most aggressive
noise reduction effect compared to the network outputs as a result of strong
regularization. However, it over-smoothed some fine structures. For example,
in Fig. 8, the vessel pointed by the green arrow was smeared out while it is
easily identifiable in NDCT as well as WGAN-VGG images. Yet, as an iterative
reconstruction method, DictRecon has its advantage over post-processing
method. As pointed by the red arrow in Fig 8, there is a bright spot which can
be seen in DictRecon and NDCT images, but is not observable in LDCT and
network processed images. Since the WGAN-VGG image is generated from LDCT
image, in which this bright spot is not easily observed, it is reasonable that
we do not see the bright spot in the images processed by neural networks. In
other words, we do not want the network to generate structure that does not
exist in the original images. In short, the proposed WGAN-VGG network is a
post-processing method and information that is lost during the FBP
reconstruction cannot easily be recovered, which is one limitation for all the
post-processing methods. On the other hand, as an iterative reconstruction
method, DictRecon algorithm generates images from raw data, which has more
information than the post-processing methods.

TABLE II: Quantitative results associated with different network outputs for Figs. 5 and 7 | Fig. 5 |  | Fig. 7
---|---|---|---
| PSNR | SSIM |  | PSNR | SSIM
LDCT | 19.7904 | 0.7496 |  | 18.4519 | 0.6471
CNN-MSE | 24.4894 | 0.7966 |  | 23.2649 | 0.7022
WGAN-MSE | 24.0637 | 0.8090 |  | 22.7255 | 0.7122
CNN-VGG | 23.2322 | 0.7926 |  | 22.0950 | 0.6972
WGAN-VGG | 23.3942 | 0.7923 |  | 22.1620 | 0.6949
WGAN | 22.0168 | 0.7745 |  | 20.9051 | 0.6759
‘1 GAN | 21.8676 | 0.7581 |  | 21.0042 | 0.6632
DictRecon | 24.2516 | 0.8148 |  | 24.0992 | 0.7631
TABLE III: statistical properties of the blue rectangle areas in Figs. 5 and 7. The values are in Hounsfield Unit (HU). | Fig. 5 |  | Fig. 7
---|---|---|---
| Mean | SD |  | Mean | SD
NDCT | 9 | 36 |  | 118 | 38
LDCT | 11 | 74 |  | 118 | 66
CNN-MSE | 12 | 18 |  | 120 | 15
WGAN-MSE | 9 | 28 |  | 115 | 25
CNN-VGG | 4 | 30 |  | 104 | 28
WGAN-VGG | 9 | 31 |  | 111 | 29
WGAN | 23 | 37 |  | 135 | 33
GAN | 8 | 35 |  | 110 | 32
DictRecon | 4 | 11 |  | 111 | 13
TABLE IV: subjective quality scores (mean$\pm$sd) for different algorithms

NDCT LDCT CNN-MSE CNN-VGG WGAN-MSE WGAN-VGG WGAN GAN DictRecon Noise
Suppression - - $4.35\pm 0.24$ $3.10\pm 0.23$ $3.55\pm 0.25$ $3.20\pm 0.25$
$2.90\pm 0.26$ $3.00\pm 0.21$ $\bm{4.65\pm 0.20}$ Artifact Reduction - -
$1.70\pm 0.28$ $2.85\pm 0.32$ $3.05\pm 0.27$ $\bm{3.45\pm 0.25}$ $2.90\pm
0.28$ $3.05\pm 0.27$ $2.05\pm 0.27$ Overall Quality $3.95\pm 0.20$ $1.35\pm
0.16$ $2.15\pm 0.25$ $3.05\pm 0.20$ $3.30\pm 0.21$ $\bm{3.70\pm 0.15}$
$3.05\pm 0.22$ $3.10\pm 0.21$ $2.05\pm 0.36$

###  III-E Quantitative Analysis

For quantitative analysis, we calculated the peak-to-noise ratio (PSNR) and
structural similarity (SSIM). The summary data are in Table II. CNN-MSE ranks
the first in terms of PSNR, while WGAN is the worst. Since PSNR is equivalent
to the per-pixel loss, it is not surprising that CNN-MSE, which was trained to
minimize MSE loss, outperformed the networks trained to minimize other
feature-based loss. It is worth noting that these quantitative results are in
decent agreement with Fig. 4, in which CNN-MSE has the smallest MSE loss and
WGAN has the largest. The reason why WGAN ranks the worst in PSNR and SSIM is
because it does not include either MSE or VGG regularization. DictRecon
achieves the best SSIM and a high PSNR. However, it has the problem of image
blurring and leads to blocky and waxy artifacts in the resultant images. This
indicates that PSNR and SSIM may not be sufficient in evaluating image
quality.

In the reviewing process, we found two papers using similar network
structures. In [38], Wolterink et al. trained three networks, i.e. GAN, CNN-
MSE, and GAN-MSE for cardiac CT denoising. Their quantitative PSNR results are
consistent with our counterpart results. And Yu et al. [39] used GAN-VGG to
handle the de-alising problem for fast CS-MRI. Their results are also
consistent with ours. Interestingly, despite the high PSNRs obtained by MSE-
based networks, the authors in the two papers all claim that GAN and VGG loss
based networks have better image quality and diagnostic information.

To gain more insight into the output images from different approaches, we
inspect the statistical properties by calculating the mean CT numbers
(Hounsfield Units) and standard deviations (SDs) of two flat regions in Figs.
5 and 7 (marked by the blue rectangles). In an ideal scenario, a noise
reduction algorithm should achieve mean and SD to the gold standard as close
as possible. In our experiments, the NDCT FBP images were used as gold
standard because they have the best image quality in this dataset. As shown in
Table III, Both CNN-MSE and DictRecon produced much smaller SDs compared to
NDCT, which indicates they over-smoothed the images and supports our visual
observation. On the contrary, WGAN produced the closest SDs yet smaller mean
values, which means it can reduce noise to the same level as NDCT but it
compromised the information content. On the other hand, the proposed WGAN-VGG
has outperformed CNN-VGG, WGAN-MSE and other selected methods in terms of mean
CT numbers, SDs, and most importantly visual impression.

In addition, we performed a blind reader study on 10 groups of images. Each
group contains the same image slice but processed by different methods. NDCT
and LDCT images are also included for reference, which are the only two
labeled images in each group. Two radiologists were asked to independently
score each image in terms of noise suppression and artifact reduction on a
five-point scale (1 = unacceptable and 5 = excellent), except for the NDCT and
LDCT images, which are the references. In addition, they were asked to give an
overall image quality score for all the images. The mean and standard
deviation values of the scores from the two radiologists were then obtained as
the final evaluation results, which are shown in Table. IV. It can be seen
that CNN-MSE and DictRecon give the best noise suppression scores while the
proposed WGAN-VGG outperforms the other methods for artifact reduction and
overall quality improvement. Also, *-VGG networks provide higher scores than
*-MSE networks in terms of artifact reduction and overall quality but lower
scores for noise suppression. This indicates that MSE loss based networks are
good at noise suppression at a loss of image details, resulting in an image
quality degradation for diagnosis. Meanwhile, the networks using WGAN give
better overall image quality than the networks using CNN, which supports the
use of WGAN for CT image denoising.

(a) VGG Map of Full Dose Image

(b) VGG Map of Quarter Dose Image

(c) Absolute Difference

Figure 9: VGG feature maps of full dose and quarter dose images in Fig. 5 and
their absolute difference.

###  III-F VGG Feature Extractor

Since VGG network is trained on natural images, it may cause concerns on how
well it performs on CT image feature extraction. Thus, we displayed two
feature maps of normal dose and quarter dose images and their absolute
difference in Fig. 9. The feature map contains 512 small images of size
$32\times 32$. We organize these small images into a $32\times 16$ array. Each
small image emphasizes a feature of the original CT image, i.e. boundaries,
edges, or whole structures. Thus, we believe VGG network can also serve a good
feature extractor for CT images.

##  IV Discussions and Conclusion

The most important motivation for this paper is to approach the gold standard
NDCT images as much as possible. As described above, the feasibility and
merits of GAN has been investigated for this purpose with the Wasserstein
distance and the VGG loss. The difference between using the MSE and VGG losses
is rather significant. Despite the fact that networks with MSE would offer
higher values for traditional figures of merit, VGG loss based networks seem
desirable for better visual image quality with more details and less
artifacts.

The experimental results have demonstrated that using WGAN helps improve image
quality and statistical properties. Comparing the images of CNN-MSE and WGAN-
MSE, we can see that the WGAN framework helped to avoid over-smoothing effect
typically suffered by MSE based image generators. Although CNN-VGG and WGAN-
VGG visually share a similar result, the quantitative analysis shows WGAN-VGG
enjoys higher PSNRs and more faithful statistical properties of denoised
images relative to those of NDCT images. However, using WGAN/GAN alone reduced
noise but at the expense of losing critical features. The resultant images do
not show a strong noise reduction. Quantitatively, the associated PSNR and
SSIM increased modestly compared to LDCT but they are much lower than what the
other networks produced. Theoretically, WGAN/GAN network is based on
generative model and may generate images that look naturally yet cause a
severe distortion for medical diagnostics. This is why an additive loss
function such as MSE and VGG loss should be added to guarantee the image
content remains the same.

It should be noted that the experimental data contain only one noise setting.
Networks should be re-trained or re-tuned for different data to adapt for
different noise properties. Especially, networks with WGAN are trying to
minimize the distance between two probability distributions. Thus, their
trained parameters have to be adjusted for new datasets. Meanwhile, since the
loss function of WGAN-VGG is a mixture of feature domain distance and the GAN
adversarial loss, they should be carefully balanced for different dataset to
reduce the amount of image content alternation.

The denoising network is a typical end-to-end operation, in which the input is
a LDCT image while the target is a NDCT image. Although we have generated
images visually similar to NDCT counterparts in the WGAN-VGG network, we
recognize that these generated images are still not as good as NDCT images.
Moreover, noise still exists in NDCT images. Thus, it is possible that VGG
network has captured these noise features and kept them in the denoised
images. This could be a common problem for all the denoising networks. How to
outperform the so-called gold standard NDCT images is an interesting open
question. Moreover, image post-denoising methods also suffer from the
information loss during the FBP reconstruction process. This phenomena is
observed in the comparison with DictRecon result. A better way to incorporate
the strong fitting capability of neural network and the data completeness of
CT data is to design a network that maps directly from raw projection to the
final CT images, which could be a next step of our work.

In conclusion, we have proposed a contemporary deep neural network that uses a
WGAN framework with perceptual loss function for LDCT image denoising. Instead
of focusing on the design of a complex network structure, we have dedicated
our effort to combine synergistic loss functions that guide the denoising
process so that the resultant denoised results are as close to the gold
standard as possible. Our experiment results with real clinical images have
shown that the proposed WGAN-VGG network can effectively solve the well-known
over-smoothing problem and generate images with reduced noise and increased
contrast for improved lesion detection. In the future, we plan to incorporate
the WGAN-VGG network with more complicated generators such as the networks
reported in [26, 27] and extend these networks for image reconstruction from
raw data by making a neural network counterpart of the FBP process.

## References

  * [1] D. J. Brenner and E. J. Hall, “Computed tomography —- an increasing source of radiation exposure,” _New England J. Med._ , vol. 357, no. 22, pp. 2277–2284, Nov. 2007.
  * [2] A. B. De Gonzalez and S. Darby, “Risk of cancer from diagnostic x-rays: estimates for the UK and 14 other countries,” _The Lancet_ , vol. 363, no. 9406, pp. 345–351, Jan. 2004.
  * [3] J. Wang, H. Lu, T. Li, and Z. Liang, “Sinogram noise reduction for low-dose CT by statistics-based nonlinear filters,” in _Med. Imag. 2005: Image Process._ , vol. 5747. International Society for Optics and Photonics, Apr. 2005, pp. 2058–2067.
  * [4] J. Wang, T. Li, H. Lu, and Z. Liang, “Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for low-dose x-ray computed tomography,” _IEEE Trans. Med. Imag._ , vol. 25, no. 10, pp. 1272–1283, Oct. 2006.
  * [5] A. Manduca, L. Yu, J. D. Trzasko, N. Khaylova, J. M. Kofler, C. M. McCollough, and J. G. Fletcher, “Projection space denoising with bilateral filtering and CT noise modeling for dose reduction in CT,” _Med. Phys._ , vol. 36, no. 11, pp. 4911–4919, Nov. 2009.
  * [6] M. Beister, D. Kolditz, and W. A. Kalender, “Iterative reconstruction methods in x-ray ct,” _Physica Medica: Eur. J. Med. Phys._ , vol. 28, no. 2, pp. 94–108, Apr. 2012.
  * [7] A. K. Hara, R. G. Paden, A. C. Silva, J. L. Kujak, H. J. Lawder, and W. Pavlicek, “Iterative reconstruction technique for reducing body radiation dose at CT: Feasibility study,” _Am. J. Roentgenol._ , vol. 193, no. 3, pp. 764–771, Sep. 2009.
  * [8] J. Ma, J. Huang, Q. Feng, H. Zhang, H. Lu, Z. Liang, and W. Chen, “Low-dose computed tomography image restoration using previous normal-dose scan,” _Med. Phys._ , vol. 38, no. 10, pp. 5713–5731, Oct. 2011.
  * [9] Y. Chen, X. Yin, L. Shi, H. Shu, L. Luo, J.-L. Coatrieux, and C. Toumoulin, “Improving abdomen tumor low-dose CT images using a fast dictionary learning based processing,” _Phys. Med. Biol._ , vol. 58, no. 16, p. 5803, Aug. 2013.
  * [10] P. F. Feruglio, C. Vinegoni, J. Gros, A. Sbarbati, and R. Weissleder, “Block matching 3D random noise filtering for absorption optical projection tomography,” _Phys. Med. Biol._ , vol. 55, no. 18, p. 5401, Sep. 2010.
  * [11] B. De Man and S. Basu, “Distance-driven projection and backprojection in three dimensions,” _Phys. Med. Biol._ , vol. 49, no. 11, p. 2463, May. 2004.
  * [12] R. M. Lewitt, “Multidimensional digital image representations using generalized Kaiser–Bessel window functions,” _J. Opt. Soc. Amer. A_ , vol. 7, no. 10, pp. 1834–1846, Oct. 1990.
  * [13] B. R. Whiting, P. Massoumzadeh, O. A. Earl, J. A. O’Sullivan, D. L. Snyder, and J. F. Williamson, “Properties of preprocessed sinogram data in x-ray computed tomography,” _Med. Phys._ , vol. 33, no. 9, pp. 3290–3303, Sep. 2006.
  * [14] I. A. Elbakri and J. A. Fessler, “Statistical image reconstruction for polyenergetic x-ray computed tomography,” _IEEE Trans. Med. Imag._ , vol. 21, no. 2, pp. 89–99, Feb. 2002.
  * [15] S. Ramani and J. A. Fessler, “A splitting-based iterative algorithm for accelerated statistical x-ray CT reconstruction,” _IEEE Trans. Med. Imag._ , vol. 31, no. 3, pp. 677–688, Mar. 2012.
  * [16] E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” _Phys. Med. Biol._ , vol. 53, no. 17, p. 4777, Aug. 2008.
  * [17] Y. Liu, J. Ma, Y. Fan, and Z. Liang, “Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction,” _Phys. Med. Biol._ , vol. 57, no. 23, p. 7923, Nov. 2012\.
  * [18] Z. Tian, X. Jia, K. Yuan, T. Pan, and S. B. Jiang, “Low-dose CT reconstruction via edge-preserving total variation regularization,” _Phys. Med. Biol._ , vol. 56, no. 18, p. 5949, Nov. 2011.
  * [19] Q. Xu, H. Yu, X. Mou, L. Zhang, J. Hsieh, and G. Wang, “Low-dose x-ray CT reconstruction via dictionary learning,” _IEEE Trans. Med. Imag._ , vol. 31, no. 9, pp. 1682–1697, Sep. 2012.
  * [20] Y. Zhang, X. Mou, G. Wang, and H. Yu, “Tensor-based dictionary learning for spectral CT reconstruction,” _IEEE Trans. Med. Imag._ , vol. 36, no. 1, pp. 142–154, Jan. 2017.
  * [21] D. Kang, P. Slomka, R. Nakazato, J. Woo, D. S. Berman, C.-C. J. Kuo, and D. Dey, “Image denoising of low-radiation dose coronary CT angiography by an adaptive block-matching 3D algorithm,” in _SPIE Med. Imag._ International Society for Optics and Photonics, Mar. 2013, pp. 86 692G–86 692G.
  * [22] G. Wang, M. Kalra, and C. G. Orton, “Machine learning will transform radiology significantly within the next 5 years,” _Med. Phys._ , vol. 44, no. 6, pp. 2041–2044, Mar. 2017.
  * [23] G. Wang, “A perspective on deep imaging,” _IEEE Access_ , vol. 4, pp. 8914–8924, Nov. 2016.
  * [24] C. Dong, C. C. Loy, K. He, and X. Tang, “Image super-resolution using deep convolutional networks,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 38, no. 2, pp. 295–307, Feb 2016.
  * [25] H. Chen, Y. Zhang, W. Zhang, P. Liao, K. Li, J. Zhou, and G. Wang, “Low-dose CT denoising with convolutional neural network,” 2016. [Online]. Available: Figures:1610.00321
  * [26] H. Chen, Y. Zhang, M. K. Kalra, F. Lin, Y. Chen, P. Liao, J. Zhou, and G. Wang, “Low-dose CT with a residual encoder-decoder convolutional neural network,” _IEEE Trans. Med. Imag._ , vol. 36, no. 12, pp. 2524–2535, Dec. 2017.
  * [27] E. Kang, J. Min, and J. C. Ye, “A deep convolutional neural network using directional wavelets for low-dose x-ray CT reconstruction,” 2016. [Online]. Available: arXiv:1610.09736
  * [28] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _Proc. Int. Conf. Med. Image Comput. Comput.-Assisted Intervention_. Springer, Nov. 2015, pp. 234–241.
  * [29] J. Johnson, A. Alahi, and L. Fei-Fei, “Perceptual losses for real-time style transfer and super-resolution,” 2016. [Online]. Available: arXiv:1603.08155
  * [30] C. Ledig, L. Theis, F. Huszar, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang, and W. Shi, “Photo-realistic single image super-resolution using a generative adversarial network,” 2016. [Online]. Available: arXiv:1609.04802
  * [31] M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein GAN,” 2017. [Online]. Available: arXiv:1701.07875
  * [32] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” in _Advances Neural Inform. Process. Syst._ , 2014, pp. 2672–2680.
  * [33] I. Goodfellow, “NIPS 2016 tutorial: Generative adversarial networks,” 2017. [Online]. Available: arXiv:1701.00160
  * [34] A. Brock, T. Lim, J. M. Ritchie, and N. Weston, “Neural photo editing with introspective adversarial networks,” 2016. [Online]. Available: arXiv:1609.07093
  * [35] J.-Y. Zhu, P. Krähenbühl, E. Shechtman, and A. A. Efros, “Generative visual manipulation on the natural image manifold,” in _Eur. Conf. Comput. Vision_. Springer, 2016, pp. 597–613.
  * [36] P. Isola, J.-Y. Zhu, T. Zhou, and A. A. Efros, “Image-to-image translation with conditional adversarial networks,” 2016. [Online]. Available: arXiv:1611.07004
  * [37] D. Nie, R. Trullo, C. Petitjean, S. Ruan, and D. Shen, “Medical image synthesis with context-aware generative adversarial networks,” 2016. [Online]. Available: arXiv:1612.05362
  * [38] J. M. Wolterink, T. Leiner, M. A. Viergever, and I. Isgum, “Generative adversarial networks for noise reduction in low-dose CT,” _IEEE Trans. Med. Imag._ , Dec. 2017.
  * [39] S. Yu, H. Dong, G. Yang, G. Slabaugh, P. L. Dragotti, X. Ye, F. Liu, S. Arridge, J. Keegan, D. Firmin _et al._ , “Deep de-aliasing for fast compressive sensing mri,” 2017. [Online]. Available: arXiv:1705.07137
  * [40] M. Arjovsky and L. Bottou, “Towards principled methods for training generative adversarial networks,” in _NIPS 2016 Workshop on Adversarial Training. In review for ICLR_ , vol. 2016, 2017.
  * [41] I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, and A. Courville, “Improved training of wasserstein GANs,” 2017. [Online]. Available: arXiv:1704.00028
  * [42] M. Nixon and A. S. Aguado, _Feature Extraction & Image Process._, 2nd ed. Academic Press, 2008.
  * [43] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” 2014. [Online]. Available: arXiv:1409.1556
  * [44] S. Srinivas, R. K. Sarvadevabhatla, K. R. Mopuri, N. Prabhu, S. S. Kruthiventi, and R. V. Babu, “A taxonomy of deep convolutional neural nets for computer vision,” _Frontiers Robot. AI_ , vol. 2, p. 36, 2016.
  * [45] AAPM, “Low dose CT grand challenge,” 2017. [Online]. Available: http://www.aapm.org/GrandChallenge/LowDoseCT/#
  * [46] D. Kingma and J. Ba, “Adam: A method for stochastic optimization,” 2014. [Online]. Available: arXiv:1412.6980
  * [47] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin _et al._ , “Tensorflow: Large-scale machine learning on heterogeneous distributed systems,” 2016. [Online]. Available: arXiv:1603.04467
  * [48] J. Deng, W. Dong, R. Socher, L. J. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in _2009 IEEE Conference on Computer Vision and Pattern Recognition_ , June 2009, pp. 248–255.


# A double-period oscillation signal in millimeter emission of the radio
galaxy NGC 1275

Pengfei Zhang  Department of Astronomy, School of Physics and Astronomy, Key
Laboratory of Astroparticle Physics of Yunnan Province, Yunnan University,
Kunming 650091, People’s Republic of China<EMAIL_ADDRESS><EMAIL_ADDRESS>Zhongxiang Wang Department of Astronomy, School of
Physics and Astronomy, Key Laboratory of Astroparticle Physics of Yunnan
Province, Yunnan University, Kunming 650091, People’s Republic of China;
<EMAIL_ADDRESS><EMAIL_ADDRESS>Shanghai Astronomical
Observatory, Chinese Academy of Sciences, 80 Nandan Road, Shanghai 200030,
China  Mark Gurwell  Center for Astrophysics $\mid$ Harvard & Smithsonian, 60
Garden Street, Cambridge, MA 02138, USA  Paul J. Wiita  Department of Physics,
The College of New Jersey, 2000 Pennington Road, Ewing, NJ 08628-0718, USA

###### Abstract

The nearby Seyfert type galaxy NGC 1275 contains a bright radio nucleus at its
center, revealed through high-spatial resolution imaging to be the source of
the jets emanating from the galaxy. Coincident with the emergence of a new
component C3 in the nucleus since 2005, flux densities from NGC 1275, at least
at radio, millimeter (mm), and $\gamma$-ray frequencies, had been increasing
up through 2017 and leveled off afterwards. We analyze the long-term light
curves of the nucleus that span the rising trend to 2015 July, and find a pair
of approximately year-long quasi-periodic oscillations, with periods of
$P_{l}\simeq 345$ d and $P_{h}\simeq 386$ d respectively, in emission at
1.3-mm wavelength. We discuss the case that there would be a long precession
period $P_{\rm prec}\simeq 9$ yr, causing the appearance of $P_{h}$ that is
slightly higher than $P_{l}$. The accretion disk around the central
supermassive black hole (SMBH) would be precessing at $P_{\rm prec}$, induced
by either the Lense-Thirring effect or the existence of a companion SMBH. In
the two scenarios, $P_{l}$ would be the jet wobbling timescale or the SMBH
binary period respectively. The finding, which could be verified through high-
spatial resolution mm imaging, would not only identify the nature of the jet
variation but also help reveal the full features of the galaxy.

##  1 Introduction

The galaxy NGC 1275, with its radio counterpart named 3C 84, is located at the
center of the Perseus cluster. As a cluster’s central galaxy that is close-by
($z=0.0176$), it has been extensively observed in various studies. The galaxy
itself has been classified as a Seyfert 1.5 type based on the emission lines
present in its optical spectrum (Véron-Cetty & Véron, 2010). In its central
active nucleus, the supermassive black hole (SMBH) has a mass of $\sim
3.4\times 10^{8}$ $M_{\odot}$, as estimated from kinematic measurements of the
surrounding molecular gas obtained from near-infrared observations (Wilman et
al., 2005). The nucleus contains jets seen at radio frequencies along the
south and north directions (Vermeulen et al., 1994). The radio core has been
further resolved to mainly consist of two components, C1 and C3, over the past
15 years (see Figure 4). While the C1 region is considered to be the base of
the jet close to the center of the galaxy, the C3 component emerged in 2005
and has been moving towards the south (Nagai et al., 2010). In addition, high-
energy and very high-energy $\gamma$-ray emission has also been detected from
NGC 1275 (Abdo et al., 2009; Aleksić et al., 2012), presumably arising from
the jets.

The core appeared to have had a rising trend in brightness at multiple
wavelengths in 2003–2017 (Hodgson et al. 2018; Britzen et al. 2019; Gulati et
al. 2021; see Figure 1). This long-term trend, as well as flare-like
variations of NGC 1275, were revealed by high-resolution imaging to primarily
arise from the emergence and flux density changes of the C3 component (Nagai
et al., 2010; Hodgson et al., 2018; Kino et al., 2018). However according to
the studies at 2–14 mm reported by Hodgson et al. (2018) for the time period
of 2013–2017, the C1 component could be brighter and more variable at 2–3 mm
than C3. Also during 2015 August–September, a “flip” in the direction of C3
was seen, followed by a monotonic flux density increase (Kino et al., 2018).
The flux density peaks in 2016–2017 seen in the long-term radio and mm light
curves (cf., Figure 1) likely were the result of the flip, thought to be
caused by the collision between the the tip of jet at C3 and a dense clump
(Kino et al., 2018).

Figure 1: OVRO 15 GHz (panel A), SMA 1.3-mm and 870-$\mu$m (black and pink
respectively; panel B), and Fermi-LAT 0.1–500 GeV (panel C) light curves. The
maximum interval of two adjacent data points in the 1.3-mm light curve is 239
d, marked by a green cross in panel B. The dashed and dash-dotted lines
indicate the times of 2015 Aug. 1 and 2016 Jan. 15. respectively (cf. Section
2.2). The radio and mm flux density peaks in 2016–2017 have been attributed to
a collision between the C3 jet and a clump of gas (Kino et al., 2018).

Related to these intriguing results for the core of the NGC 1275 jets, hints
of precession of the jets have been reported. The derived precession
timescales, however, have a very wide range: from as long as $\sim 33$ Myr to
explain the observed multi-kpc scale radio (and X-ray) bubbles (Dunn et al.,
2006); to $\sim 40$–100 yr from parsec-scale morphological changes and long-
term flux density variations of the jets (Britzen et al., 2019); $\sim 28.8$
yr, from determining the projected position angle of the C1 component (Dominik
et al., 2021); or $\sim$6 yr, from a possible helical path of C3 (with respect
to C1; Hiura et al. 2018).

Since the multiwavelength emission from the core of the jets in NGC 1275 has
been bright, monitoring of it has made long-term light curve data available at
several wavelengths. We took advantage of these data, and conducted searches
for possible quasi-periodic oscillation (QPO) signals for a further study of
the jets. Surprisingly, a double-period oscillation signal at 1.3 mm
wavelength was found. Previously, a search for QPO signals was conducted with
the 9-yr $\gamma$-ray data obtained using the Fermi Gamma-ray Space Telescope
(Fermi), but no periodicities were found (Nemmen et al., 2020).

In the following Section 2, we describe the data and periodicity search
analysis and provide the periodicity results. In Section 3, the origin of the
periodicities is discussed. Additional information about the source and
results of our analyzes are provided in Appendices A and B.

##  2 Data analysis and Results

###  2.1 Data description and Fermi $\gamma$-ray data analysis

The archival data we used to search for periodicity signals include 1.3-mm and
870-$\mu$m monitoring of NGC 1275 with the Submillimeter Array (SMA; Gurwell
et al. 2007), 15 GHz monitoring with the Owens Valley Radio Observatory (OVRO;
Readhead et al. 1989; Richards et al. 2011), and 0.1–500 GeV monitoring with
the Large Area Telescope (LAT; Atwood et al. 2009) onboard Fermi.

The SMA has been monitoring over 400 sources on a several-day cadence at these
two wavelengths since 2002, and NGC 1275 (3C 84) is one of them. The two flux
density light curves obtained span over 18 yr, approximately from 2002 Nov. to
2021 Jan./Feb.. For the 1.3-mm light curve (Figure 3), the median and mean
time intervals between two adjacent data points are 3.13 and 9.35 d
respectively, and the maximum interval is $\sim$239 d around MJD 52808 (2003
June 18). On many days two or three measurements were made. For the 870-$\mu$m
light curve (Figure 1), the sampling cadence is much sparser compared to that
of the 1.3-mm one: the median, mean, and maximum time intervals between two
adjacent data points are 15.47, 60.24, and 876.65 d, respectively.

The OVRO 40-m telescope has been monitoring over 1800 AGN on a $\sim$3 d
cadence since 2008. NGC 1275 (J0319+4130 in the OVRO target list) is one of
the targets. The flux density light curve (Figure 1) we used is approximately
10 yr long, from 2008 Jan. 8 (MJD 54473) to 2017 Nov. 15 (MJD 58072). The
median and mean intervals between two adjacent data points are 3.71 and 5.36
d, respectively, and the maximum interval is 84.75 d.

Figure 2:  Upper left panel: WWZ power for the 1.3-mm light curve in
2002–2015-July. Two periodicities, around $\sim 345$ d and $\sim 386$ d, are
clearly seen. Upper right panel: LSP power (red line) and summed WWZ power
(black line) obtained for the 1.3-mm light curve. The best fit to the
underlying power spectral density is shown as the blue dash-dotted line. The
5$\sigma$ significance curve, obtained from the light curve simulation, is
shown as the green dashed curve. The power peaks for the two periodicities are
above the 5$\sigma$ significance curve. After taking into account a trial
number of 350, the 345 d (386 d) periodicity still has a significance greater
than (close to) 5$\sigma$. Bottom panel: Comparison of LSP powers from the
1.3-mm light curve in the time periods from 2002 to 2021 (black dashed), 2002
to 2016 Jan. 15 (green), and 2002 to 2015 Aug. 01 (yellow). These power peaks
are nearly the same, though the highest one for the 346 d periodicity is from
the 2002 to 2016 Jan. 15 time period.

For the _Fermi_ -LAT data, we selected the 0.1-500 GeV Pass 8 _Front+Back_
SOURCE class photon-like events from 2008 Aug. 4 to 2021 Apr. 16 in a region
of $20^{\circ}\times 20^{\circ}$ centered at the position of NGC 1275. In the
fourth _Fermi_ Large Area Telescope catalog (4FGL; Abdollahi et al. 2020), the
source is named 4FGL J0319.8+4130 with the coordinate R.A.=$03^{h}19^{m}49\fs
8$, Decl.=$41^{\circ}30^{{}^{\prime}}43\farcs 57$. The events were reduced by
selecting the zenith angle
$\rm<90^{\circ},~{}DATA\\_QUAL>0,and~{}LAT\\_CONFIG=1$ to obtain high-quality
data in the good time intervals. We built a model file that contained the
spectral parameters of all known 4FGL sources in the selected region. The
spectral form of NGC 1275 in the 4FGL is log-parabola, described by the form
$dN/dE=K(E/E_{0})^{-\alpha_{s}-\beta_{s}\rm{log}(E/E_{b})}$. The binned
maximum likelihood analysis was first applied to the whole data, and for NGC
1275, we obtained an average photon flux of $\rm(341.7\pm 3.1)\times
10^{-9}~{}photons\,cm^{-2}\,s^{-1}$ in the 0.1–500 GeV band with an overall
test statistic (TS) value of $\sim$147500\. The corresponding best-fit
spectral parameters we obtained are $\alpha_{s}=2.078\pm 0.005$,
$\beta_{s}=0.054\pm 0.003$, and $E_{b}=0.953\pm 0.010$ GeV, in good agreement
with that reported in 4FGL.

The results from the above analysis were saved as a new model file, based on
which the maximum likelihood analysis was performed again for constructing a
0.1–500 GeV light curve of NGC 1275 (Figure 1). We chose 15 d as the time-bin
size, so that the TS values obtained from each time-bin data are mostly $\geq
25$.

###  2.2 Periodicity search and determination

Two methods, the Weighted Wavelet Z-transform (WWZ; Foster 1996) and the
generalized Lomb-Scargle Periodogram (LSP; Lomb 1976; Scargle 1982;
Zechmeister & Kürster 2009) were employed to search for periodicities in the
radio, mm, and $\gamma$-ray light curves. Only in the 1.3-mm light curve did
we detect a significant QPO signal, which actually consists of two
periodicities (upper panels of Figure 2). This QPO signal is present
throughout the entire light curve (i.e., from 2002 to 2021). However, since
the broad flux density peaks in 2016–2017 at radio and mm bands were likely
due to the jet-clump collision in 2015 Aug.–Sept., we also analyzed only the
2002–2015 Aug. 01 light curve data, and the resulting power peaks are nearly
the same as those from the whole data (bottom panel of Figure 2). In addition,
we found that the power peak at the lower period value actually strongest when
the data restricted to the slightly longer period of 2002–2016 Jan. 15 are
used. In any case, because of the evidence for a “flip” seen in high-
resolution imaging (Kino et al., 2018), in the remainder of this paper we
focus on the results from the 2002–2015 July data.

The period values of the two periodicities were obtained by fitting the two
LSP power peaks with Gaussian functions. The values are $P_{l}=344.6\pm 10.0$
d and $P_{h}=386.0\pm 10.9$ d, where the uncertainties were taken as the half-
widths at half-maximum of each peak.

In other light curves, weak signals were seen in the 870-$\mu$m light curve,
but they did not have significances sufficiently high to make any QPO claims
(see Appendix B and Figure 5). We note that the 870-$\mu$m light curve only
has 111 data points, and is much sparser than the 1.3-mm light curve, so the
lack of a significant signal is not surprising.

Figure 3: Top panel: Long-term 1.3-mm light curve of NGC 1275. The fits given
in Section 2.2 were used to detrend the long-term variations in the light
curve, which are a degree 3 polynomial (red and blue dashed lines) and a
constant (horizontal blue dashed line) respectively for the rising portion
from 2002 to 2017 and the part after the peak from 2017 to 2021. Bottom panel:
Detrended light curve, obtained from dividing the observed light curve by the
best-fits to the underlying trends shown in the top panel. The double-
sinusoidal best-fit model given by the LSP code is shown as the yellow dashed
curve.

We evaluated the significances of the two periodicities by simulating the
light curve. We used a smoothly bending power law plus a constant,
$P(f)=Af^{-\alpha}[1+(f/f_{bend})^{\beta-\alpha}]^{-1}+C$, to model the
underlying power spectrum density (PSD) of the flux density variations
(González-Martín & Vaughan, 2012). The model function consists of the
parameters $A$, $\alpha$, $\beta$, $f_{bend}$, and $C$, which are the
normalization, low frequency slope, high frequency slope, bend frequency, and
Gaussian noise, respectively. A maximum likelihood method (Barret & Vaughan,
2012) was used to fit the PSD, and we obtained the best-fit values $A=0.27\pm
0.17$, $\alpha=0.42\pm 0.32$, $\beta=0.59\pm 0.29$, $f_{bend}=0.012\pm 0.109$,
and $C=0.78\pm 0.47$, respectively, so for this source a single power-law
would be acceptable. Using the model, we generated $3\times 10^{7}$ light
curves (Emmanoulopoulos et al., 2013) and obtained the PSD of each.
Significance curves were built by counting the PSD data points at each
frequency considered. The 5$\sigma$ significance curve obtained in this
fashion is shown in Figure 2. The two power peaks of the periodicities reach
above the significance curve. After considering the trial number in our
temporal analysis, which is 350, $P_{l}$ and $P_{h}$ have only modestly
reduced significances, which respectively are greater than and close to
5$\sigma$.

The LSP code provides the best-fit model for the QPO, which consists of two
sinusoidal functions in a form of $A\sin[2\pi f(t-t_{0})]+A_{0}$. The
parameters for $P_{l}$ ($P_{h}$) are $A=1.47\pm 0.13$ ($1.09\pm 0.13$),
$f=2.900\pm 0.016\times 10^{-3}$ d-1 ($2.585\pm 0.020\times 10^{-3}$ d-1),
$t_{0}=137.51\pm 4.68$ ($48.89\pm 7.45$), $A_{0}=5.37\pm 0.09$ ($5.21\pm
0.09$). In order to clearly show the QPO variations in the 1.3-mm light curve,
we detrended its long-term flux density changes. The rising part before 2015
August 01 (of duration MJD 52577–57235) was fitted with a degree 3 polynomial,
$2.4-4.4\times 10^{-5}t+5.2\times 10^{-7}t^{2}-1.4\times 10^{-11}t^{3}$ Jy
(with $t$ in days after MJD 52577). The possible QPO variations, obtained by
dividing the observed light curve with this function, are shown in Figure 3.
The normalized sinusoidal functions are overplotted. The model can
approximately describe the variations. In order to understand if the
variations in the 2016–2021 part of the light curve could be related, we
extended this model fitting to the end of the light curve. The part
corresponding to the 2016–2017 variation peak (MJD 57235–57870) was still
quite well fit by the polynomial, and because the remaining part of 2017–2021
(MJD 57870–59258) is relatively flat, we fit it with a constant and normalized
by the constant (9.18 Jy). The data points in general do not show a clear
agreement with the double-sinusoidal modulation. The comparison may explain
the analysis result that the inclusion of the 2016–2021 light curve part do
not help increase the significance.

As the SMA observations are prohibited from pointing to targets close to the
Sun and thus the light curve we analyzed has a certain pattern that might
induce artifactual signals (VanderPlas, 2018), we conducted detailed analyses
to confirm the real presence of the QPO signal. The analyses are described in
Appendix section C. We did not find any evidence that the signal was an
artifact.

##  3 Discussion

The QPO signal we have found is the first time one has been seen at mm
wavelengths (to our knowledge) among the QPOs reported in emission at
different wavelengths from Active Galactic Nuclei (AGN). While the radio and
$\gamma$-ray light curves show similar long-term variation patterns, we did
not detect significant signals in them, in particular, the 15 GHz one (Figure
1). In this case, the mm observations may be able to more sensitively catch
the variations of the jet’s emission (Hodgson et al., 2018), thus revealing
the QPO signal.

Also uniquely, the QPO signal consists of two close periodicities and thus
presents a very intriguing case. The thus-far reported AGN QPOs are widely
discussed to reflect instances such as periodic variations of the accretion
disk surrounding a central SMBH or the orbital period of a binary SMBH (see,
e.g., Gierliński et al. 2008; King et al. 2013; Ackermann et al. 2015; Zhang &
Wang 2021), and often compared to relatively well-studied cases in stellar-
mass black hole systems as well. In the latter, QPOs are seen to appear as a
series of harmonically related peaks (Ingram & Motta, 2019), $\nu_{0}$,
2$\nu_{0}$, …, (n+1)$\nu_{0}$, where $\nu_{0}$ is the fundamental peak
frequency. The double periods in our case obviously do not fit in this type of
QPO phenomena, as $P_{l}$ and $P_{h}$ correspond to frequencies of 3.36$\times
10^{-8}$ Hz and 3.00$\times 10^{-8}$ Hz respectively. Theoretical simulations
have shown that a thick disk in an AGN surrounding a SMBH may undergo acoustic
p-mode oscillations when excited by an external perturbation (Rubio-Herrera &
Lee, 2005). Such an oscillation signal could consist of multiple frequencies,
which may have been seen in a couple of periodic variability cases found in
AGN (e.g., An et al. 2013). However, the dominant frequencies should follow
the ratio of 2:3, which is not the case here.

The disk/jet precession scenario due to the Lense-Thirring effect (or frame-
dragging effect) has been applied to the explanations of QPO phenomena in
X-ray binaries (e.g., Stella & Vietri 1998; Ingram & Motta 2019), tidal
disruption events (e.g., Pasham et al. 2019), and blazars (e.g., Bhatta &
Dhital 2020). In this scenario, a thick disk around a rapidly spinning SMBH
could be tilted and such a tilted disk would precess at a timescale of $\sim
0.02(M/10^{8}M_{\sun})(r/r_{g})^{3}$ d (Merritt et al., 2010; Bhatta & Dhital,
2020), where $M$ and $r_{g}=GM/c^{2}$ are the mass and the gravitational
radius, respectively, of the SMBH. For the NGC 1275 case, $M=3.4\times 10^{8}\
M_{\sun}$ (Wilman et al., 2005), and if the inner part of its presumed tilted
disk is at $\sim 20\ r_{g}$, the timescale would be $\sim 1$ yr, possibly
leading to similarly precessing jets (Liska et al., 2018). So the QPO signal
we have found to arise from the jets of NGC 1275 could be another case of this
scenario. However, there are two close periodicities in our case. We note that
in recent general relativistic magnetohydrodynamic simulations (Liska et al.,
2018), the jets were found to be wobbling around the disk at an approximate
timescale of $10^{4}r_{g}/c$, which is $\sim 0.5$ yr for the case of NGC 1275.
If we consider that $P_{l}$ arises from the jets and is due to wobbling, the
disk precession timescale, $P_{\rm prec}$, could be much longer (and stable),
resulting in the weaker periodicity $P_{h}$ through a beat frequency where
$P_{\rm prec}=1/(1/344.6\,{\rm d}-1/386.0\,{\rm d})\simeq 9$ yr. This $P_{\rm
prec}$ value would put the inner part of the disk in NGC 1275 at $\sim
40r_{g}$. This scenario, with a jet wobbling at $P_{l}$ and precessing at
$P_{\rm prec}$, is thus a possible way to explain the observed two
periodicities. Whether such a wobbling jet could induce a significant QPO
signal is not clear and should be further investigated.

In addition, we note the well-studied case of the stellar-mass black hole
binary system SS 433, which has an orbital period of $\sim 13$ d and whose
jets are precessing with a period of $\sim 162$ d (Margon, 1984). Quite a few
X-ray binaries, including SS 433, show a period much longer than the known
orbital period, and this is considered to possibly reflect a precessing disk
that surrounds the central compact star in such an X-ray binary (Wijers &
Pringle, 1999). So we could consider $P_{l}$ as an orbital period of a binary
SMBH system and take the jets to be precessing at $P_{\rm prec}$ estimated
above. The ratio of $P_{\rm prec}$ to $P_{l}$ would be $\simeq 9.3$, similar
to that (162 d/13 d$\simeq 12.5$) in SS 433. Such ratios of other X-ray
binaries are generally in a range of 10–100 (Wijers & Pringle, 1999). This
similarity leads to the possible scenario where the galaxy contains jets
precessing with a $\sim 9$ yr period and that the center could consist of two
SMBHs orbiting each other at a period of $\sim 345$ d.

There is evidence that NGC 1275 experienced a merger with another galaxy
probably $\sim$300 Myr ago (Carlson et al., 1998; Conselice et al., 2001). A
binary SMBH system would be a natural result of the merger. However current
studies show that it often takes giga-year for the two SMBHs from each galaxy
to form a compact binary after a merger event (e.g., Yu 2002). Considering 345
d as the orbital period, the binary would have a short separation distance of
$\sim$0.003 pc (here assuming equal masses for the two SMBHs, with
$M_{1}=M_{2}$ and $M_{1}+M_{2}=3.4\times 10^{8}\,M_{\odot}$). Given the
compactness (Yu, 2002), either this putative binary SMBH might not be the
result of the previous merger but may have formed from an earlier merger, or
the formation of a compact SMBH binary could simply be a factor of a few
faster than theoretically expected.

As shown by Hodgson et al. (2018), the C1 component was brighter and more
variable at 2–3 mm wavelengths than the C3 component before the apparent 2015
August C3 collision with a clump of gas. Thus the QPO signal we have found
could very well arise from C1. Indeed, by analyzing the position angle changes
of C1 at 43 GHz, Dominik et al. (2021) found a $\sim 29$ yr precession or a
$\sim 21$ yr precession plus a $\sim 3.4$ yr nutation. Although these possible
periods do not match those of the QPO, it should be noted that their
observations were at radio and over a limited time span of 9 yr. As mm
observations may be more attuned to catching the variations of the jets in NGC
1275, long-term mm high-resolution imaging of the jets would potentially be
able to determine the precession period (or periods), verify our QPO results,
and confidently identify the nature of the jet variation.

We thank the anonymous referee for helpful suggestions, Y. Luo for discussions
about binary supermassive black hole evolution, and F. Xie for discussions
about QPOs in X-ray binaries. The Submillimeter Array is a joint project
between the Smithsonian Astrophysical Observatory and the Academia Sinica
Institute of Astronomy and Astrophysics and is funded by the Smithsonian
Institution and the Academia Sinica. This research has made use of data from
the OVRO 40-m monitoring program which was supported in part by NASA grants
NNX08AW31G, NNX11A043G and NNX14AQ89G, and NSF grants AST-0808050 and
AST-1109911, and private funding from Caltech and the MPIfR. This research is
supported by the National Key R & D Program of China under grant No.
2018YFA0404204, the joint foundation of Department of Science and Technology
of Yunnan Province and Yunnan University $[$2018FY001 (-003)$]$, and the
National Natural Science Foundation of China (12163006, 11633007). Z.W.
acknowledges the support by the Original Innovation Program of the Chinese
Academy of Sciences (E085021002).

## References

  * Abdo et al. (2009) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJ, 699, 31
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33
  * Ackermann et al. (2015) Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 813, L41
  * Aleksić et al. (2012) Aleksić, J., Alvarez, E. A., Antonelli, L. A., et al. 2012, A&A, 539, L2
  * An et al. (2013) An, T., Baan, W. A., Wang, J.-Y., Wang, Y., & Hong, X.-Y. 2013, MNRAS, 434, 3487
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071
  * Barret & Vaughan (2012) Barret, D., & Vaughan, S. 2012, ApJ, 746, 131
  * Bhatta & Dhital (2020) Bhatta, G., & Dhital, N. 2020, ApJ, 891, 120
  * Britzen et al. (2019) Britzen, S., Fendt, C., Zajaček, M., et al. 2019, Galaxies, 7, 72
  * Carlson et al. (1998) Carlson, M. N., Holtzman, J. A., Watson, A. M., et al. 1998, AJ, 115, 1778
  * Conselice et al. (2001) Conselice, C. J., Gallagher, John S., I., & Wyse, R. F. G. 2001, AJ, 122, 2281
  * Dominik et al. (2021) Dominik, R. M., Linhoff, L., Elsässer, D., & Rhode, W. 2021, MNRAS, 503, 5448
  * Dunn et al. (2006) Dunn, R. J. H., Fabian, A. C., & Sanders, J. S. 2006, MNRAS, 366, 758
  * Emmanoulopoulos et al. (2013) Emmanoulopoulos, D., McHardy, I. M., & Papadakis, I. E. 2013, MNRAS, 433, 907
  * Foster (1996) Foster, G. 1996, AJ, 112, 1709
  * Gierliński et al. (2008) Gierliński, M., Middleton, M., Ward, M., & Done, C. 2008, Nature, 455, 369
  * González-Martín & Vaughan (2012) González-Martín, O., & Vaughan, S. 2012, A&A, 544, A80
  * Gulati et al. (2021) Gulati, S., Bhattacharya, D., Bhattacharyya, S., et al. 2021, MNRAS, 503, 446
  * Gurwell et al. (2007) Gurwell, M. A., Peck, A. B., Hostler, S. R., Darrah, M. R., & Katz, C. A. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 375, From Z-Machines to ALMA: (Sub)Millimeter Spectroscopy of Galaxies, ed. A. J. Baker, J. Glenn, A. I. Harris, J. G. Mangum, & M. S. Yun, 234
  * Hiura et al. (2018) Hiura, K., Nagai, H., Kino, M., et al. 2018, PASJ, 70, 83
  * Hodgson et al. (2018) Hodgson, J. A., Rani, B., Lee, S.-S., et al. 2018, MNRAS, 475, 368
  * Ingram & Motta (2019) Ingram, A. R., & Motta, S. E. 2019, New A Rev., 85, 101524
  * King et al. (2013) King, O. G., Hovatta, T., Max-Moerbeck, W., et al. 2013, MNRAS, 436, L114
  * Kino et al. (2018) Kino, M., Wajima, K., Kawakatu, N., et al. 2018, ApJ, 864, 118
  * Liska et al. (2018) Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2018, MNRAS, 474, L81
  * Lister et al. (2018) Lister, M. L., Aller, M. F., Aller, H. D., et al. 2018, ApJS, 234, 12
  * Lomb (1976) Lomb, N. R. 1976, Ap&SS, 39, 447
  * Margon (1984) Margon, B. 1984, ARA&A, 22, 507
  * Merritt et al. (2010) Merritt, D., Alexander, T., Mikkola, S., & Will, C. M. 2010, Phys. Rev. D, 81, 062002
  * Nagai et al. (2010) Nagai, H., Suzuki, K., Asada, K., et al. 2010, PASJ, 62, L11
  * Nemmen et al. (2020) Nemmen, R., de Menezes, R., & Paschalidis, V. 2020, in Perseus in Sicily: From Black Hole to Cluster Outskirts, ed. K. Asada, E. de Gouveia Dal Pino, M. Giroletti, H. Nagai, & R. Nemmen, Vol. 342, 167–171
  * Pasham et al. (2019) Pasham, D. R., Remillard, R. A., Fragile, P. C., et al. 2019, Science, 363, 531
  * Readhead et al. (1989) Readhead, A. C. S., Lawrence, C. R., Myers, S. T., et al. 1989, ApJ, 346, 566
  * Richards et al. (2011) Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29
  * Rubio-Herrera & Lee (2005) Rubio-Herrera, E., & Lee, W. H. 2005, MNRAS, 357, L31
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835
  * Stella & Vietri (1998) Stella, L., & Vietri, M. 1998, ApJ, 492, L59
  * Suzuki et al. (2012) Suzuki, K., Nagai, H., Kino, M., et al. 2012, ApJ, 746, 140
  * VanderPlas (2018) VanderPlas, J. T. 2018, ApJS, 236, 16
  * Vermeulen et al. (1994) Vermeulen, R. C., Readhead, A. C. S., & Backer, D. C. 1994, ApJ, 430, L41
  * Véron-Cetty & Véron (2010) Véron-Cetty, M. P., & Véron, P. 2010, A&A, 518, A10
  * Wijers & Pringle (1999) Wijers, R. A. M. J., & Pringle, J. E. 1999, MNRAS, 308, 207
  * Wilman et al. (2005) Wilman, R. J., Edge, A. C., & Johnstone, R. M. 2005, MNRAS, 359, 755
  * Yu (2002) Yu, Q. 2002, MNRAS, 331, 935
  * Zechmeister & Kürster (2009) Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577
  * Zhang & Wang (2021) Zhang, P., & Wang, Z. 2021, ApJ, 914, 1

##  Appendix A Illustration of the jets in the nucleus of NGC 1275

In Figure 4, the 15 GHz radio image, taken in 2014 from the Monitoring Of Jets
in Active galactic nuclei with VLBA Experiments (MOJAVE; Lister et al. 2018)
is shown. In the right panel of the figure, the C1 and C3 components are
indicated. Previously there was a relatively faint component C2, which was
located west (at the right side) of C3 (Suzuki et al., 2012).

Figure 4: Left: 15 GHz radio image of the central 20 mas region of NGC 1275.
Right: enlarged view of the central part, with the C1 and C3 components
indicated.

##  Appendix B Long-term light curves of NGC 1275 and periodicity analysis
results

In Figure 1, the long-term radio, mm, and $\gamma$-ray light curves of NGC
1275 are shown. We employed the WWZ and LSP analysis methods to search for
periodicities in the radio, 870-$\mu$m, and $\gamma$-ray light curves, and the
results are shown in Figure 5. The time periods were from the beginning of the
data to 2015 August 01, before the probable C3-clump collision. No significant
signals in the three bands were found. We also conducted the analysis of the
entire light curve data for each of the three datasets, and the results were
the same.

Figure 5: WWZ powers and LSP powers for the 15 GHz radio (top left),
870-$\mu$m (top right), and Fermi-LAT 0.1–500 GeV (bottom) light curves of NGC
1275. For each of these bands there is one modest power peak close to the
periodicity ranges determined from the 1.3-mm light curve (shown with pink and
blue shadings). In addition, the 870 $\mu$m band shows a rather strong peak
around 500 d and the others have some signals longer than 500 d. However,
considering the sparseness of the data points at 870 $\mu$m and high values of
each underlying power spectral density at large periods, none of the signals
at 870-$\mu$m or $\gamma$-rays are significant.

##  Appendix C Analysis checks for the QPO signal in NGC 1275

As shown in Figure 6, the 1.3-mm light curve of NGC 1275 has an unavoidable
gap in May of each year when the source is within a 25 degree radius from the
Sun (the solar avoidance zone for SMA observations), which might cause
artifact signals in analyses for detecting periodic signals. In order to
investigate whether this observational pattern could induce the signal we have
found, we conducted the following analyses as additional checks for the
signal. Here we include the light curve data after 2015 August 01, as the gaps
are present throughout the entirity of the data.

Firstly we conducted the same analysis to the 1.3-mm flux density light curves
of J0359+509 (R. A.=$03^{h}59^{m}29\fs 7472$,
Decl.=$50^{\circ}57^{{}^{\prime}}50\farcs 161$) and J2253+161 (R.
A.=$22^{h}53^{m}57\fs 7479$, Decl.=$16^{\circ}08^{{}^{\prime}}53\farcs 560$),
obtained with the SMA in the same monitoring program. The former is a source
close to our target in the sky and the latter was one of the most frequently
observed with the SMA, and so provides a particularly good test for production
of any artifactual signal. For J0359+509, a paucity of data points in
approximately April–May of each year is present probably for the same reason
that NGC 1275 could not be observed then, but in total there are 287 data
points in its light curve (Figure 7). For J2253+161, there is an obvious gap
in March of each year but its light curve has 919 data points (Figure 8). From
the WWZ and LSP analyses, no significant signals that were similar to those
appearing in the power spectra of NGC 1275 were seen.

We further investigated how the gap might affect the QPO signal in NGC 1275.
Examining the yearly light curve (of the detrended one), we noted that there
were only a few data points in 2002–2004 and 2021. Ignoring these sparse data
points, the LSP power peak of $P_{l}$ would be lowered from $\sim 50$ to $\sim
40$, while the power peak of $P_{h}$ would be kept nearly the same (Figure 6).
We estimated the mean number of the data points per month ($\simeq$4) based on
the 2005–2020 light curve, and obtained its standard deviation ($\sigma_{d}$).
We then generated randomized pseudo-data points to fill in the May gaps in
this light curve. The times of these data points were randomly chosen in each
May and their fluxes were randomly chosen within the $\pm 2\sigma_{d}$ range
(the average flux uncertainty was given to each generated flux). The generated
data points (64 in total) are shown in the top right panel of Figure 6. Using
the two functions for detrending the long-term variations, we brought back the
generated data points to fit into the original light curve for 2005–2020. We
then conducted a LSP analysis of this gap-filled light curve and the resulting
power spectrum density (PSD) is shown in the bottom panel of Figure 6. Both
the $P_{l}$ and $P_{h}$ periodicity peaks were lowered to $\sim 30$. We
modeled the PSD and then generated 5$\times 10^{5}$ light curves based on the
model fit. The 3$\sigma$ and 4$\sigma$ significance curves were obtained from
the simulated light curves. The two power peaks remain significant at a
$>3\sigma$ level.

We understand the above results as follows. As the May gaps are the major
pattern seen in the observed light curve, we tested to remove the pattern by
filling the generated data points. The LSP analysis results still show that
the two-periodicities are significant, although the significances are lowered,
particularly the $P_{l}$ one. Since we generated approximately 9% of the
observed data points, for which their fluxes were randomly chosen from the
$\pm 2\sigma_{d}$ range around the smoothly fitted underlying trends, we
effectively added noise to any putative QPO signal present in the light curve.
Thus the significances naturally became lower than those in the observed light
curve, but the signal remains strong. We note that a weak peak in the LSP PDS
at $\sim 57$ d is slightly increased instead, now reaching nearly the
4$\sigma$ significance curve (Figure 6). The reason for this change is not
clear, but in any case we mostly modified the low-frequency structures of the
PSD by adding the data points and the high-frequency structures might not be
significantly affected.

From these analyses, by searching for periodic signals from a nearby source
and a most frequently observed source whose light curve has a similar gap
pattern, and by filling the gaps in the light curve of NGC 1275 with random
data points, we conclude that the two periodicities we have found are real
signals, not artifacts induced by the constrained observation patterns.

Figure 6: Top left: detrended 1.3-mm light curve of NGC 1275. A gap in May of
each year (marked by green dotted lines) is clearly seen. In addition, there
are only a few data points in 2002–2004 and 2021 (shown in red). Top right:
gap-filled light curves for NGC 1275 in 2005–2020. For each gap in 2005-2020,
we randomly generated data points based on the mean number of data points per
month (i.e., setting times in each gap) and $\pm 2\sigma_{d}$ range (i.e.,
setting fluxes at those times). The generated data points are shown in yellow
in the yearly folded light curve (bottom panel). Adding back the values
calculated from the detrending functions to the generated data points, we
obtained the gap-filled 2005–2020 light curve. Bottom: power curves from LSP
analyses of the original 1.3-mm light curve (yellow), 2005–2020 light curve
(orange), gap-filled 2005–2020 light curve (black). The best fit to the black
power curve is shown as the blue dash-dotted curve, based on which the
3$\sigma$ (green dashed) and 4$\sigma$ (yellow dotted) significance curves
were obtained.

Figure 7: Top figures: 1.3-mm SMA light curves of J0359+509, which is close to
NGC 1275 in the sky, where the left one shows the standard view of the
observations and the right one shows the yearly light curve, illustrating a
gap in approximately April–May of each year (similar to that of NGC 1275).
Bottom figure: WWZ power (left panel) and LSP power and summed WWZ power (red
and black lines respectively in the right panel) obtained for the 1.3-mm light
curve. No signals of QPOs, real or artifactual, are present.

Figure 8: Top figures: 1.3-mm SMA light curves of J2253+161, a very well
measured source, where the left one shows the standard view of the
observations and the right one shows the yearly light curve. There is a gap in
March of each year. Bottom figure: WWZ power (left panel) and LSP power and
summed WWZ power (red and black lines respectively in the right panel) for the
1.3-mm light curve. No signals could be identified with a sufficiently high
significance and, critically, none are present at the periods seen in the
light curve of NGC 1275.