\section{Introduction and related work}
\label{sec:intro}

Noise is inevitable in biological and medical imaging. Equipment imperfection in microscopy leads to image artifacts, some signal frequencies are undersampled in magnetic resonance scans, exposure time is restricted to capture dynamic information, or to avoid sample bleaching and phototoxicity. Therefore, image analysis in medicine and biology often relies on image reconstruction. The reconstruction's primary objective is to approximate "clean" image signal $s = x + s'$ by recovering absent information $s'$ from the acquired image $x$, or decomposing the acquisition $x = s + n$ into signal $s$ and degrading noise $n$. Given that the ways to approach both reconstruction tasks are similar, we will further use the terms "reconstruction", "restoration" and "denoising" interchangeably.

Deep learning advances have improved image restoration~\cite{lehtinen2018noise2noise, krull2019noise2void, krull2019probabilistic, buchholz2019cryo, batson2019noise2self, weigert2018content}. Convolutional neural networks are currently outperforming traditional (non-trainable) methods such as non-local means~\cite{buades2011non} or BM3D~\cite{dabov2007image} in open benchmarks~\cite{lemarchand2019opendenoising}, \emph{e.g}\onedot} \def\Eg{\emph{E.g}\onedot~in fluorescent microscopy denoising~\cite{zhang2019poisson}. From the deep learning point of view, the most intuitive way to approach the task is to train a neural network to reconstruct clean target from the noisy image. This approach is often referred to as Noise2Clean. During the training, network parameters are optimized to minimize the error between predicted signal $\hat{s}$ and true signal $s$. But acquiring clean ground truth is challenging, often requiring multiple measurements and their subsequent registration ~\cite{weigert2018content, zhang2019poisson}.

Perhaps surprisingly, training is also possible without clean ground truth. One way to train a denoising neural network in such way is to use two independently corrupted versions (i.e. multiple acquisitions) of the same image, as was proposed by Lehtinen~\emph{et al}\onedot~\cite{lehtinen2018noise2noise}. Noise2Noise reuses corrupted versions of the same image both as input and target to effectively reconstruct the clean signal across various noise types and image domains including magnetic resonance imaging (MRI) scans~\cite{lehtinen2018noise2noise}. This approach is mathematically justified by an assumption that the network learns the expected value of target pixels (mean, median or mode — depending on used loss function). One of its main limitations is the requirement to have two independently corrupted copies of the data.
While it is relatively straightforward to synthetically corrupt ground truth images twice, acquiring two independent measurements is usually not feasible.

To mitigate this requirement, multiple methods were developed to train a denoising network from a single corrupted measurement. Noise2Void by Krull~\emph{et al}\onedot~\cite{krull2019noise2void} and Noise2Self by Batson and Royer~\cite{batson2019noise2self} exploited the idea that signal is not statistically independent across the image, but noise is, so self-supervision could be applied. Noise2Void is trained to predict a central pixel under the blind-spot from a patch around it, while Noise2Self generalizes the approach by introducing an additional masking procedure. Laine~\emph{et al}\onedot~\cite{laine2019improved} further improved blind-spot networks and performed on par with Noise2Noise.
Noisier2Noise proposed by Moran~\emph{et al}\onedot~\cite{moran2019noisier2noise} outperformed Noise2Void in self-supervised denoising but required a predefined model of the noise. Finally, Noise2Void has been extended into a probabilistic model by Krull~\emph{et al}\onedot~\cite{krull2019probabilistic} by estimating the parameters of the noise distribution from either camera characteristics measurement or noise-free copy of the data, which allowed to compete with supervised models on a publicly available microscopy dataset~\cite{zhang2019poisson}.

There are opportunities to improve these models by using recent discoveries about neural network architectures, loss functions and training procedures. However, while extended approaches described above dominated traditional training-free denoising algorithms, they could not outperform supervised methods like Noise2Noise which use more information. This motivates investigating ways to use more data during training instead. In health imaging and cell biology, it is common to have a stack of 2D planes that could be combined in a single volumetric image. For MRI, these planes are resonance signals from a certain sample depth~\cite{lustig2008compressed}, for cell microscopy, acquisitions from consecutive focal planes~\cite{dalgarno2010multiplane}. Given that each image consists of multiple planes, it may also be possible to take advantage of mutual information between neighboring planes and denoise an image stack from self (see Figure~\ref{fig:ssim-show}). This approach was applied to denoise cryogenic electron microscopy (cryo-EM)~\cite{buchholz2019cryo}, low-dose computer tomography (CT) and MRI data~\cite{wu2019consensus}. There, a stack of tomograms was split into even and odd planes, and a network was trained to predict even planes using odd and vice versa. This approach is beneficial when neither clean data nor the second noisy copy is available.
On the other hand, the number of training examples in this case decreases to the number of available image stacks, which often is a strict limiting factor in neural network training.

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.45\linewidth]{fig/plane_residual_ssim_50.jpg}
\hspace{10pt}
\includegraphics[width=0.45\linewidth]{fig/plane_residual_ssim_80.jpg}

\end{center}
   \caption{Structural similarity (SSIM) of neighboring planes in MRI dataset. Planes 50 and 51 from one stack, followed by the residual between neighboring planes; and the same for planes 80 and 81. All images are scaled in range $[-0.5, 0.5]$. Residuals have near-zero means and standard deviations 0.05 and 0.07 respectively. SSIM between planes 50 and 80 is $0.38$. High structural similarity suggests that we can use plane neighbors for reconstruction, but within limited distance.}
\label{fig:ssim-show}
\end{figure}

In this work, we present Noise2Stack, an extension of the Noise2Noise model that makes use of information from the neighbouring planes to improve the quality of the image denoising (see Figures~\ref{fig:title},~\ref{fig:microscopy} for example results).
Our experiments demonstrate that Noise2Stack trained on neighbouring planes in addition to the copy of the target image, consistently outperforms the original Noise2Noise method. Moreover, we show that even without noisy copy of the target image, its performance matches or even exceeds the one of Noise2Noise method. Such training scheme allows learning from a single stack of 2D images without requiring the second independently corrupted copy, which is often hard or impossible to acquire.

\section{Methods}

\subsection{Noise2Stack}
\label{sec:n2s}

Our first main contribution is the Noise2Stack method. It is a sampling algorithm for training denoising neural networks which selects a set of planes from a provided stack on each step. For example, in order to denoise the $i$-th plane in the stack, the network can be given information from the neighbors $i+1$ and $i-1$. Below we describe the nuances of our algorithm's work.

Noise2Stack can operate in two modes. In the  copy-supervised mode, a copy of the reconstructed plane is passed as an input along with its neighbors. Without adding neighbouring planes, this is exactly equivalent to the original Noise2Noise strategy.
In the second, self-supervised mode, the reconstructed plane is held out as a training target and not used as an input. For even numbers of input planes, the self-supervised mode was used. For example, with two-plane training, to denoise the $i$-th plane, only its neighbors $i+1$ and $i-1$ were used. For odd numbers of planes, the copy-supervised mode was used: the $i$-th plane was also a part of the input.

In both modes, we train a single denoising neural network for all the planes ignoring their absolute spatial location in a stack as in Noise2Noise. Note that for the self-supervised mode, we do not require two independently corrupted copies of a stack. Instead, we use just a single noisy stack, predicting each image plane using its neighbors. For the marginal planes, we replace their absent neighbors with copies of their nearest neighbor plane (\emph{e.g}\onedot} \def\Eg{\emph{E.g}\onedot~for the first plane we use two copies of the second plane). Figure~\ref{fig:model} illustrates the difference between Noise2Noise and Noise2Stack in a self-supervised mode.

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.85\linewidth]{fig/models_wide.jpg}
\end{center}
   \caption{Noise2Noise and Noise2Stack training strategies. In Noise2Noise, two independently corrupted copies of each image plane are used (blue copy as an input, red copy as a target). In Noise2Stack we can choose whether we show reconstructed plane to the network as an input along with its neighbors (two neighbors training illustrated) or not. It gives an opportunity to train the network when the second corrupted copy of the stack is not available. Unused image planes are excluded from the training process at this step.}
\label{fig:model}
\end{figure}

\subsubsection{Neural network architecture}
\label{sec:nn}

Denoising can be done with different network architectures which often are domain-specific.
The previously proposed Noise2Noise~\cite{lehtinen2018noise2noise} and Noise2Void~\cite{krull2019noise2void} approaches are not architectures \textit{per se}, but training strategies, which can be applied with any suitable neural network backbone. Noise2Stack is no exception to this rule.
The only architectural aspect one may need to modify when adopting a new backbone is the number of input channels, which in general case is equal to the number of planes sampled from the stack at each training step.
In our experiments, we used previously published variations of a U-Net~\cite{ronneberger2015u} that were applied to similar tasks.
For MRI denoising, we used an architecture presented in the Noise2Noise paper~\cite{lehtinen2018noise2noise}, and for microscopy denoising we adopted an architecture from CSBDeep framework~\cite{weigert2018content}. Both of them were implemented using PyTorch deep learning library~\cite{paszke2019pytorch} in Python, and both are described in Supplementary material (Tables~\ref{tab:supp-network-microscopy},~\ref{tab:supp-network-mri}).

\subsubsection{Training}
\label{sec:training}

In all the experiments the neural network was trained for 100 epochs with batch size 16 on a single NVIDIA V100 GPU. As in~\cite{lehtinen2018noise2noise}, we did not use any regularization techniques, such as dropout or batch normalization. The network was trained using Adam~\cite{kingma2017adam} optimizer with initial learning rate $10^{-3}$ and default parameters $\beta_1=0.9$, $\beta_2=0.999$, $\epsilon = 10^{-8}$ without weight decay. The learning rate was decreasing every epoch according to the cosine schedule~\cite{loshchilov2016sgdr}. We used mean squared error (MSE) loss function as recommended in~\cite{lehtinen2018noise2noise}. The loss was computed directly on the network output without additional post-processing steps.

\subsection{Microscopy dataset}
\label{sec:dataset}

To test our approach on multiplane microscopy images, we acquired a novel dataset~\cite{kenny_roberts_2020_4114086}. CHO, U2OS, and RPE1 cell lines were grown in standard medium, seeded on plates, fixed with formaldehyde, stained with a fluorescent dye (Hoechst33342) and imaged with PerkinElmer Phenix high-throughput confocal microscope in fluorescence and brightfield modalities using a 20x objective. For the fluorescent modality, images were acquired in low-exposure (20 ms) and high-exposure (100 ms) modes. For the brightfield, we used a single exposure (100 ms). Higher exposure time usually translates into better image quality, alleviating Poisson noise, however, it leads to sample degradation, and lower measurement speed.

In total, three wells on a plate were imaged, one well per cell line, each with nine fields of view. Each field of view was imaged across ten focal planes separated by 2 micrometers in z-stack, resulting in 270 microscopy images of size $2160 \times 2160$ pixels. We used the first two wells (CHO and U2OS, 180 images) as a training set, the first three fields of view from the last well (RPE1, 30 images) as a validation set, and the last six fields of view (RPE1, 60 images) as a test set.

\subsection{Noise generation and post-processing for MRI}
\label{sec:noise-gen}
\label{sec:post}

To generate noise for MRI dataset, we followed the procedure presented in~\cite{lehtinen2018noise2noise}. We used a Bernoulli process to sample from the frequency space of the image signal by performing a Fourier transform, and retaining each frequency $k$ with probability $p(k) = e^{-\lambda |k|}$. Signal from selected frequencies was weighted with the inverse of their respective probabilities, and omitted frequencies were replaced with zeros. As in~\cite{lehtinen2018noise2noise}, the sampling parameter $\lambda$ was chosen to preserve 10\% of frequencies relative to a full sampling. Resulting spectra were then transformed back to the original space by inverse Fourier transform.

The noise coming from spectrum undersampling provides an opportunity for post-processing.
For MRI images, we transformed the network's input and output into frequency space, and copied all the  sampled frequencies from the input image spectrum to the output as is. After copying frequencies, we transformed the resulting combined spectrum back to the image with inverse fast Fourier transform. As was shown in~\cite{lehtinen2018noise2noise}, enforcing the exact preservation of input frequencies in such a way improves the result of image restoration. This post-processing was done under the assumption that we need to restore only those frequencies which were masked. As such post-transformation is possible only when the noise comes from frequency undersampling and is not feasible for \emph{e.g}\onedot} \def\Eg{\emph{E.g}\onedot~Poisson or Gaussian noise in microscopy, we report the performance both before and after post-processing for MRI.

\subsection{Metrics}
\label{sec:metrics}

We used peak signal-to-noise ratio (PSNR) as a primary metric to compare our results with previous reports~\cite{lehtinen2018noise2noise}. In addition, we measured structural similarity index (SSIM) and normalized root mean squared error (NRMSE). For fluorescent microscopy images, before measuring metrics we applied percentile normalization and affine transformation as described in~\cite{weigert2018content}. First, we percentile-normalized ground truth between $0.1$ and $99.9$, then we affinely scaled and translated predictions to minimize MSE between predictions and ground truth. Data range for metric calculation was set to 1.

\section{Experiments}
\label{sec:res}

We evaluated our Noise2Stack method against several baselines on two multiplane datasets from different biomedical domains. First, we tested the method's ability to remove the synthetic noise from MRI in the full IXI brain scan T1 dataset.\footnote{http://biomedic.doc.ic.ac.uk/brain-development/downloads/IXI/IXI-T1.tar} Noise was generated by undersampling frequencies as described in \ref{sec:noise-gen}. Second, we denoised multiplane low-exposure microscopy images in brightfield and fluorescent modalities. For the fluorescent images, we evaluated the denoising results against registered high-exposure images and in the downstream segmentation task.

\subsection{MRI}
\label{sec:res-mri}

 MRI dataset consists of volumetric images of $256 \times 256$ pixels with 150 planes each, from 60 subjects in total. We selected 100 middle planes from the full stack since marginal planes contain less information and are not as informative for noise measurements. We selected 4800 planes from 48 subjects for the training set, 200 planes from 2 subjects for the validation set and 1000 planes from 10 subjects for the test set. The list of subjects selected for testing matches the validation set from~\cite{lehtinen2018noise2noise}. Our validation set is a small subset of the training set from~\cite{lehtinen2018noise2noise} and is used to select the best model based on the validation loss prior to testing. During pre-processing pixel intensities from each image were scaled and shifted to range $[-0.5, 0.5]$, no additional normalization was applied. As an only data augmentation technique, we randomly translated images in the range $[0, 64]$ pixels .

\label{sec:res-baseline}
We measured the baseline performance on the whole test set for direct inverse fast Fourier transform reconstruction of undersampled spectra as in~\cite{lehtinen2018noise2noise}. Our test set had an average PSNR of $20.7$ dB for directly reconstructed images. In addition, we directly applied the post-processing step by combining the spectra of input and output images. For the copy-supervised Noise2Noise training, this is possible because we assume that we have two independently corrupted copies of data. We consider this noise level more appropriate to evaluate improvement in image quality for strategies that use two image copies. When two input images spectra are combined, mean PSNR rises to $22.4$ dB (Table~\ref{tab:results}).

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.9\linewidth]{fig/prediction_crop.jpg}
\end{center}
   \caption{Results of Noise2Stack MRI denoising compared to Noise2Noise, input image and reference ground truth. Peak signal-to-noise ratio in bottom right corner is calculated against reference image. }
\label{fig:title}
\end{figure}

\label{sec:res-n2n}
To ensure the correctness of the following experiments and architecture implementation, we replicated the results of Noise2Noise paper~\cite{lehtinen2018noise2noise}. Here, independently corrupted copies of a single image plane were used as an input and as an output. A single plane Noise2Noise model increases PSNR to $29.1$, and to $30.9$ with further post-processing. Our test set results for the original Noise2Noise training strategy are slightly ($0.8$ dB) worse than previously reported on the validation dataset~\cite{lehtinen2018noise2noise}.

\label{sec:res-n2s}
Adding four neighbors to the Noise2Noise input increased PSNR by $1.8$ dB both before and after post-processing (Table~\ref{tab:results}). Overall, limiting the number of neighbouring planes to four (two from each side) proved optimal, resulting in the greatest performance benefits with reasonable computational costs in both modes (Figure~\ref{fig:psnr-ssim}). Using a single image stack for training in a self-supervised mode with four neighbours in the input, we outperform the original Noise2Noise model by $0.8$ dB.

Next, we investigated whether position in the stack affected denoising performance. We evaluated the best trained model across different plane locations and found that the marginal planes have higher PSNR after the reconstruction as well as for the baseline (Figure~\ref{fig:to-baseline}A).

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.99\linewidth]{fig/metric_vs_planes_used.jpg}
\end{center}
  \caption{Noise suppression metrics dependent on the number of planes used in training for magnetic resonance imaging (MRI) data. Peak signal-to-noise ratio averaged across the test set images (PSNR~$\uparrow$; panel A), structural similarity (SSIM~$\uparrow$; panel B), normalized root mean squared error (NRMSE~$\downarrow$; panel C; all y-axis) for increasing number of planes used in the training (x-axis) using a copy-supervised (blue), or self-supervised (orange) approach with post-processing (solid line) and without (dashed line), training details are disclosed in text.
   Vertical bars: 95\% confidence intervals. }
\label{fig:psnr-ssim}
\end{figure}

\begin{table}[h]
\centering
\caption{Trained model performance on magnetic resonance imaging (MRI) and fluorescence cell microscopy datasets. For MRI, metrics were measured against clean image copy without synthetic noise. For microscopy, metrics were measured against high-exposure sample which we considered clean. For each model average results before and after post-processing are shown alongside.
}
\label{tab:results}
\begin{tabular}{lllcccccc} \toprule
    Dataset & Denoising method  & $N$ input & \multicolumn{2}{c}{PSNR $\uparrow$}  & \multicolumn{2}{c}{SSIM $\uparrow$} & \multicolumn{2}{c}{NRMSE $\downarrow$} \\
    && planes & before & after & before & after & before & after \\ \midrule
      MRI & - & -            & 20.7   & 22.4  & 0.38   & 0.43  & 0.244  & 0.201 \\ \cmidrule{2-9}

    & Noise2Noise &     1      & 29.1   & 30.9  & 0.85   & 0.87  & 0.092  & 0.075 \\ \cmidrule{2-9}
    & Noise2Stack &    1+2            & 30.3   & 32.2  & 0.88   & 0.90  & 0.080   & 0.064 \\
    & (copy-supervised) &    1+4            & 30.9   & 32.7  & 0.89   & 0.91  & 0.075  & 0.061 \\
    &  &    1+6            & 30.8   & 32.6  & 0.89   & 0.91  & 0.075  & 0.061 \\ \cmidrule{2-9}
    &  Noise2Stack&    2            & 28.6   & 31.3  & 0.84   & 0.87  & 0.097  & 0.072 \\
    &  (self-supervised)&    4            & 29.2   & 31.7  & 0.85   & 0.88  & 0.091  & 0.069 \\
    &  &    6            & 29.0     & 31.5  & 0.85   & 0.87  & 0.093  & 0.070 \\ \midrule
    Microscopy & - & - & 25.8 && 0.40 && 0.388 & \\ \cmidrule{2-9}
     & BM3D & 1 & 31.0 && 0.68 && 0.222 & \\ \cmidrule{2-9}
     & Noise2Clean & 1 & 32.8 && 0.77 && 0.174 & \\ \cmidrule{2-9}
    & Noise2Void & 1 & 26.8 && 0.44 && 0.357 & \\ \cmidrule{2-9}
    & Noise2Stack & 2 & 31.6 && 0.73 && 0.195 & \\
    & (self-supervised) & 4 & 31.6 && 0.73 && 0.197 & \\ \bottomrule
\end{tabular}

\end{table}

\subsection{Microscopy}

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.4\linewidth]{fig/metric_vs_plane_location.jpg}
\hspace{10pt}
\includegraphics[width=0.4\linewidth]{fig/sanger_metric_vs_plane_location_w_bm3d.jpg}
\end{center}
   \caption{The impact of plane location on model performance.
   \textbf{A:} Average peak signal-to-noise ratio on the test set of MRI dataset (y-axis), ranging the target plane location across the 100 available ones (x-axis). Red: baseline (direct inverse Fourier transform); blue: copy-supervised Noise2Stack with five input planes (two neighbors from each side); error band: 95\% confidence interval for 10 subjects.
   \textbf{B:} Average peak signal-to-noise ratio on the test set of fluorescent microscopy dataset (y-axis), ranging the target plane location across the 10 available ones (x-axis). Red: baseline (low-exposure image); other colors: models in comparison; error band: 95\% confidence interval for 6 stacks.}
\label{fig:to-baseline}
\end{figure}

Microscopy dataset consists of 270 pairs of low-exposure and high-exposure fluorescent modality images and 270 respective brightfield images. We split them 180/30/60 between training/validation/test sets. The detailed dataset description provided in Section~\ref{sec:dataset}. Before the training, camera background pattern was calculated as a median across all images in each modality, and subtracted from each image before normalization between $3$ and $99.8$ percentiles as in~\cite{weigert2018content}. Background subtraction was done to remove correlated signal and prevent the network from learning identity mapping (Figure~\ref{fig:supp-background}). During training, images were randomly rotated by \{0, 90, 180, 270\} degrees, flipped both horizontally and vertically, and cropped to $256 \times 256$ pixels.

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.98\linewidth]{fig/microscopy_predictions_show_s4_p2_w_bm3d.jpg}
\end{center}
   \caption{Microscopy denoising results compared to low-exposure and high-exposure images in fluorescent and brightfield modalities (the second plane from the fourth test stack was taken as an illustrative example). Segmentation of fluorescent images with a pre-trained neural network illustrates denoising effects on downstream analysis; red: probability map thresholded by $0.5$ and cleaned.}
\label{fig:microscopy}
\end{figure}

We compared Noise2Stack with supervised denoising method, Noise2Clean, which uses low-exposure images as inputs and high-exposure images as training targets. The second stack of low-exposure images was not available, therefore, training Noise2Noise network was not possible. In addition, we evaluated our strategy against another self-supervised method, Noise2Void~\cite{krull2019noise2void}, which does not require a noise model, does not put any restrictions on volumetric data properties, and treats all the planes independently of their position in the stack. For the training, we used the default Noise2Void configuration recommended by the authors in their GitHub repository. As a baseline, we used the unsupervised BM3D denoising algorithm~\cite{dabov2007image, makinen2019exact}.

High-exposure fluorescent images were used as a ground truth when denoising low-exposure images in supervised way and calculating metrics. Brightfield images naturally exhibit higher levels of noise comparing to fluorescence, thus, 100 ms exposure images could not possibly be used as clean targets. Therefore, we only evaluated brightfield denoising results visually.

Figure~\ref{fig:to-baseline}B illustrates peak signal-to-noise ratio for different models as well as for raw images with respect to the plane position in a stack. Experiment results are summarized in Table~\ref{tab:results}. Unlike MRI data, four-plane input did not show considerable improvement comparing to two-plane mode. Thus, in both Figures~\ref{fig:to-baseline} and~\ref{fig:microscopy} we showed denoising results obtained using two neighbouring planes as input.
Figure~\ref{fig:microscopy} shows visual comparison of the denoising results with low-exposure (noisy) and high-exposure (clean) microscopy images. Noise2Void has shown only moderate improvement comparing to the original images adding $1$ dB to the baseline PSNR of $25.8$ dB. At the same time Noise2Stack was able to capture the signal in a self-supervised way, achieving PSNR $31.6$ dB and producing visually almost as good results as supervised Noise2Clean, although $1.2$ dB less by PSNR. Both Noise2Stack and Noise2Clean outperformed the BM3D baseline in our experiments.

Denoising is usually done during pre-processing to ease the downstream image analysis for both machine algorithms and human researchers. To demonstrate the applicability of Noise2Stack as a general first step in the image analysis pipeline, we segmented fluorescence microscopy images with a pre-trained neural network. For pre-training we used 3072 high-exposure fluorescent images from another experiment with cell nuclei stained by the same protocol and annotated in a semi-automated way~\cite{fishman2019segmenting}. From each image we subtracted its mean value, then divided by the standard deviation. Images were segmented independently of their position in the stack one by one without any aggregation, because small nuclei could appear in lower planes and vanish in upper ones or vice versa. We did not produce manual annotation for our test data and thus did not quantify segmentation results. Instead, we visually evaluated the consistency of predictions and the amount of artifacts introduced by noise. In Figure~\ref{fig:microscopy} we show semantic segmentation probability maps for low-exposure, denoised, and high-exposure images. Segmentation of low-exposure images denoised with Noise2Clean and Noise2Stack produces plausible results, sometimes even better than segmentation of high-exposure images. We observed the biggest impact of denoising algorithms in marginal planes, where signal-to-noise ratio was extremely low. More examples from other test images in different stack positions are included in the Supplementary material (Figure~\ref{fig:supp-microscopy}).

\section{Discussion}

The Noise2Stack method was inspired by the observation of neighboring plane similarity in a number of volumetric biomedical datasets (Figure~\ref{fig:ssim-show}). For example, for MRI brain scan dataset that we used in our experiments (Section~\ref{sec:dataset}), the average structural similarity between neighboring planes is $0.86$. We showed that U-Net is able to capture this mutual information and effectively use it for denoising. In both copy-supervised (where original Noise2Noise belongs) and self-supervised modes, two neighbors from each side (four total) were optimal for denoising of MRI. Two neighbors, one from each side, were best for the fluorescent microscopy. Adding more neighbors to the stack did not affect the performance in copy-supervised mode and made it slightly worse for self-supervised mode, where only neighbors were used as the network's input. We hypothesize that this effect depends on the physical distance between imaged planes, and that the optimal number of planes differs between datasets. However, it is clear that further planes provide less useful information for the restoration and could lead to network overfitting.

We observed the marginal planes in MRI data to have higher PSNR for both baseline and denoised images (Figure~\ref{fig:to-baseline}). Exactly the opposite is true for fluorescent microscopy images: middle planes in a stack have the highest peak signal-to-noise ratio.
And although these patterns differ across the domains, Noise2Stack has shown stable improvement of the performance comparing to methods that do not use information from the neighbouring planes.

Unsupervised BM3D algorithm showed good results for the middle microscopy planes, however failed in the marginal ones with lower signal-to-noise ratio (Figure~\ref{fig:to-baseline}B). In complex conditions, it produced a large number of artifacts that corrupted the downstream segmentation results (Figure~\ref{fig:microscopy}). In our experiments, it was also 700 times slower than all of the deep learning methods due to the large image size and lack of GPU speed up.

Visually, the level of detail in MRI denoised with Noise2Stack in copy-supervised mode has noticeably increased comparing to Noise2Noise, as can be seen in Figure~\ref{fig:title} (additional illustrations provided in Supplementary material, Figure~\ref{fig:supp-showcase}). Brain gyri in the reconstructed images are clean and sharp even for the cerebellum where they are tiny. In self-supervised mode Noise2Stack results are expectedly less detailed and match Noise2Noise by sight. However, this mode demands twice less data, being able to learn from a single copy. For the fluorescent microscopy, Noise2Stack in self-supervised mode did not outperform supervised Noise2Clean, but has shown comparable results.

\section{Conclusion}
\label{sec:conclusion}

We introduced Noise2Stack — an approach to improve Noise2Noise results for volumetric data. We use neighboring image planes as an additional source of information for image restoration, and demonstrate that this boosts performance. If the target plane is not included, denoising is a result of the middle plane prediction. Further, we show that by using only neighboring planes it is possible to achieve Noise2Noise level of performance with a single image stack.
As the network is trained to predict the expected value of each pixel,
the method can be used even if there is no corrupted stack copy available, reducing imaging times, or avoiding additional assumptions about the noise characteristics.
Noise2Stack can be applied in all domains where measurements are made in a spatially structured manner. The limitations of sampling to secure the required amount of shared signal between acquisitions is a subject of further studies.

Our second major contribution is a dataset of matched low exposure time (noisy) and high exposure time (clean) fluorescence cell image stacks from different cell lines accompanied by the brightfield images. These freely available data, together with the results presented here, can serve as the benchmark for future methods development in the field.

\begin{ack}

This work was funded by PerkinElmer Cellular Technologies (VLTAT19682), Wellcome Trust (206194), Estonian Research Council (IUT34-4), and Estonian Centre of Excellence in IT (EXCITE) (TK148). We thank High Performance Computing Center of the Institute of Computer Science at the University of Tartu for the provided computing power.

\end{ack}

{\small
\bibliographystyle{ieee_fullname}


\section{Introduction}

Active galactic nuclei (AGN) in extreme minima or maxima states provide us with
important insight into the physics of the black hole and accretion disk region.
The group of narrow-line Seyfert 1 galaxies (NLS1 hereafter) is known to exhibit
particularly strong flux and spectral variability in X-rays (reviews
by \citet{Gallo2018a} and \citet{Komossa2018}, and references therein).
Mrk 335 \citep{Markarian1977} is a member of the class of NLS1
galaxies at redshift $z=0.0258$ with a black hole mass of $M_{\rm BH} = 2.7 \times 10^7$ M$_{\odot}$
\citep{Grier2012}. Its optical spectrum is characterized by strong high-ionization coronal
emission lines with ionization states of up to [FeX] \citep{Grupe2008} that
imply the presence of a strong EUV-X-ray continuum.
Mrk 335 was detected as a bright and variable
UV emitter in early International Ultraviolet Explorer (IUE) and Hubble
Space Telescope (HST) observations \citep[e.g.,][]{Dunn2006, Longinotti2013}.

For a time, this was one of the X-ray brightest AGN
\citep{Tananbaum1978, Pounds1987, Grupe2001, Bianchi2001}. Then the X-ray flux of Mrk 335 dropped dramatically
in 2007 \citep{Grupe2007, Grupe2008}.
Multiple epochs of exceptional flaring and dipping on its way to and beyond this
low state motivated several deeper spectroscopic observations with XMM-Newton, Suzaku, NuSTAR, and the HST
\citep{Grupe2012, Gallo2013, Kara2013, Longinotti2013, Parker2014, Gallo2015,
Wilkins2015, WilkinsGallo2015, Komossa2017, Gallo2018, Gallo2019,
Longinotti2019, Choudhury2019, Parker2019}.
These observations have shown that Mrk 335 is well fit with varying contributions of blurred reflection
and (ionized) absorption.
The reflection model explains diminished X-ray continuum emission with changes in a corona
that has collapsed in toward the black hole, and sometimes forms a collimated outflow in
X-ray flare states.
The absorption model explains several epochs of flares and fades with changes in ionization
state and covering factor of the absorbers, which were detected with the XMM-Newton
reflection grating spectrometer (RGS) and in the UV with the HST.

We report the emergence of Mrk 335 from its recent ultra-deep low state
(which lasted from late 2017 to at least January 2020), based on 45
new observations with the Neil Gehrels Swift observatory
\citep{Gehrels2004} each in all seven wavebands from optical to X-rays. These
are part of our long-term monitoring of Mrk 335, which started in May 2007. When Mrk 335 was observed on 2020 May 16, it was found already at a flux level that exceeded all observations since February 2018.
In particular, it was very bright in the UV as well and appeared at a similar UV flux
level as seen last during the flare at the end of 2011. Because of this outburst in
May-June 2020, we started an intensive monitoring campaign with Swift, and we report the results here.
We use a cosmology \citep{Wright2006} with
$H_{\rm 0}$=70 km\,s$^{-1}$\,Mpc$^{-1}$, $\Omega_{\rm M}$=0.3 and $\Omega_{\rm \Lambda}$=0.7 throughout this paper.

\section{Data analysis}

\subsection{Swift XRT}
We have monitored Mrk 335 with Swift for more than 13 years since May 2007 with a cadence of typically 1-14 days, depending on the brightness and variability state of Mrk 335. Long gaps of several months occur each year when Mrk 335 is unobservable with Swift because of its proximity to the Sun. Mrk 335 came out of the Swift Sun constraint on 16 May 2020 \citep{Grupe2020} and we resumed our monitoring (Fig. \ref{fig:lc-Swift} and Fig. \ref{fig:lc-zoom}).
The Swift X-ray telescope \citep{Burrows2005} was operating in photon-counting mode \citep{hill2004}
with typical exposure times of 1-2 ks (Tab. \ref{tab:obs-log}).

X-ray count rates were determined using the online XRT product tool at the Swift data centre in Leicester \citep{Evans2007}{\footnote{ \url{https://www.swift.ac.uk/user_objects/}}}.
For a dedicated spectral analysis of select epochs, source photons were extracted within a circle of radius 59\arcsec , and
the background was determined in a nearby circle with a radius of 236\arcsec.
During the faintest states, spectral fitting of individual spectra is not possible, and we instead measured hardness ratios, defined as $HR = \frac{hard-soft}{hard+soft}$, where {\it soft} and {\it hard} are the counts in the 0.3-1.0 and 1.0-10 keV energy bands, respectively. $HR$ was determined by applying the Bayesian estimation of hardness ratios, BEHR \citep{park06}{\footnote{ \url{http://hea-www.harvard.edu/AstroStat/BEHR/}}}.

In order to obtain a high-state spectrum for spectral fitting of a better signal-to-noise ratio (S/N), we combined the 2020 June 23 -- July 1 Swift data when Mrk 335 was brightest (Target ID 13544, segments 011-016).
For comparison, we created a low-state spectrum when Mrk 335 was in its lowest state with data from May 2019 (Target ID 33420, segments 218-221).
These are the two primary data sets used for spectral fitting.
In both cases we created new ancillary response files (arfs) by adding the arfs of the single spectra weighted by their exposure time using the {\sc{ftool}} command {\it addarf}. The coadded X-ray spectra in the band (0.3-10) keV were then analyzed with the software package {\sc xspec} \citep[version 12.10.1f;][]{Arnaud1996}.

\subsection{Swift UVOT}
Mrk 335 was observed with the UV-optical telescope \citep[UVOT;][]{Roming2005}
either in W2 or in all six filters
in order to obtain spectral energy distribution (SED) information of this rapidly varying AGN.
Observations in each filter were coadded using the task {\em{uvotimsum}}.
Source counts in all filters were then extracted in a circular region with a radius
of 5\arcsec ~, and the background was selected in a nearby region of radius 20\arcsec.
The background-corrected counts were then converted into VEGA magnitudes and
fluxes based on the latest calibration as described in \citet{Poole2008} and \citet{Breeveld2010}.
The task {\em{ uvotsource}} was used to measure the magnitudes and flux densities.
A Galactic reddening correction was applied to the
UVOT data, with a value of $E_{\rm{(B-V)}}$=0.035 \citep{Schlegel1998}, with a correction factor in each filter
according to Eq. (2) of \citet{Roming2009} and making use of the reddening curves of \citet{Cardelli1989}.

\begin{figure}

\includegraphics[clip, trim=1.3cm 5.6cm 1.3cm 2.6cm, width=\columnwidth]{fig1}
    \caption{Long-term Swift XRT and UVOT light curves of Mrk 335 from 2007 May
    to 2020 July.
    From top to bottom: XRT count rate, hardness ratio, and UVOT Vega magnitudes corrected for
    Galactic extinction. The vertical red lines mark the beginning of each second year from
    2008 to 2020.}
\label{fig:lc-Swift}
\end{figure}

\section{Results}

\subsection{Light curves and DCF analysis}
Fig.\ref{fig:lc-Swift} displays the long-term light curve of Mrk 335 in X-rays as
well as in the optical/UV since the start of our Swift monitoring program in May
2007, when it was discovered in an unusually low X-ray flux state \citep{Grupe2007}.
 Figure \ref{fig:lc-zoom} highlights the 2020 evolution in all bands.

Mrk 335 exhibited several short-duration X-ray flares over the last 13 years
when it was in a low state overall. The new rise in X-ray flux
after the long period of very low activity of more than two years is remarkable; it is
unprecedented in Mrk 335. The systematic rise
in the UV flux since about mid-2018 is even more remarkable.

Fig. {\ref{fig:lc-zoom} shows the Swift results since we resumed our monitoring program on 2020 May 16.
This light curve shows multiple rapid flare events in X-rays, including one rapid
rise by a factor of $\sim$5 within a week, and another rise by a
factor of 2 within a day. Changes in the X-ray spectral shape as shown
in the hardness ratio light curve are apparent. The XRT count
rate before May 2020 was at a very low level around 0.05 cts s$^{-1}$ for more than
two years (and it is still variable, therefore not host dominated),
and now (at the peak of the flare in late June 2020) it reached a value that is a factor of 10 higher.
The brightest state in the UV (W2) was detected on 2020 May 31 with a Galactic-reddening
corrected magnitude of 12.68\pl0.04. This is the highest state since the UV flare seen in December 2011.

In order to search for correlations and time delays between different bands during the 2020 outbursts,
a discrete correlation function \citep[DCF;][]{Edelson1988} analysis of the UVW2 and B and X-ray
light curves was carried out following \citet{Gallo2018} with Monte Carlo confidence intervals.
The UV band is closely correlated with the optical band, such that the optical is following
and the UV leading
by $\Delta t \simeq 1.5$d with a confidence interval of 3 to 0 days at the 99.99\% level.
However, no correlation between the UV and X-rays is found.

\begin{figure}
\includegraphics[trim=1.5cm 5.6cm 1.3cm 2.6cm, width=9.1cm]{fig2}
    \caption{Swift outburst light curve of Mrk 335 since 2020 May 16 (units as in Fig. \ref{fig:lc-Swift}).
The horizontal dashed red line marks an X-ray count rate of 0.05 cts/s (the lowest recorded during the 2018-2019
low state was 0.03 cts/s), while the horizontal dashed blue line marks the lowest value of the W2 magnitude
measured with Swift (in 2018). The vertical red lines mark 2020 June 1 and July 1.
    }
    \label{fig:lc-zoom}
\end{figure}

\subsection{X-ray spectral analysis}

Spectra were first fit with single power laws of photon index $\Gamma_{\rm X}$
adding Galactic absorption with a column density $N_{\rm H, Gal}=3.96\times10^{20}$ cm$^{-2}$ \citep{kalberla05}.
Because this model did not provide good spectral fits and left strong systematic residuals in all fits (Fig. \ref{fig:xray_spec}), we then applied either an ionized or neutral partial covering absorber (Tab. \,\ref{tab:spec-fits}).
The single-component partial covering model consistently requires a lower column density and covering fraction at  high state. Alternatively, the high-state spectrum is well fit by an ionized absorber.
The 2019 low-state spectrum can also be well described by a simple powerlaw plus black
body model, where the black body is taken to be representative of a soft excess, which
is dominated by a number of emission lines detected with the RGS in
low state \citep{Longinotti2019, Parker2019} and not individually resolved with Swift.
We have added such a component (with its parameters fixed to the low-state values)
when fitting the high-state spectrum, but it does not provide a significant change in
the overall fit results.
Based on the simple partial covering absorber fit to the 2020 high state, a (0.3-10) keV X-ray
luminosity of $8.2 \times 10^{43}$ erg/s is obtained, or an X-ray Eddington ratio
of $L_{\rm X}/L_{\rm Edd} \approx 0.02$, while for the 2019 low state, the intrinsic X-ray
luminosity is a factor of 1.6 lower when the same model description is adopted.

Spectral fits were also performed on a 2020 low state and intermediate state (Fig. 3). However, as
Mrk 335 varies strongly during the 2020 time interval (unlike during the deep May 2019 low state),
spectral fits are less reliable and still come with large errors, and at any given time
interval a mix of cold and ionized absorbers may contribute.
These cannot be disentangled with Swift, however.
Applying the simple partial covering model
to the 2020 low state,
we find  $N_{\rm H} = 11.32^{+4.62}_{-2.97}$ cm$^{-2}$
, while for the intermediate state, we find $N_{\rm H} = 4.32^{+1.22}_{-0.88}$ cm$^{-2}$.

\subsection{Hardness ratio variability}
We inspected the long-term (since 2007) and short-term (2020) $HR$ variability of
Mrk 335 (Fig. \ref{fig:HR}) including all data points until 2020 July 19.
Three regimes are apparent: a softer-when-fainter pattern at lowest count rate $CR$,
a softer-when-brigther trend at intermediate $CR$, and an almost constant $HR$ at highest $CR$.
These three regimes have also been identified in other absorbed AGN that are intrinsically bright \citep[e.g.,][]{Connolly2014}.
Results of a regression analysis (available in {\sc{midas}}) of the three variability regimes are shown in Fig. \ref{fig:HR}. A correlation analysis at $CR <0.1$ cts/s gives $r_s = +0.5, T_S = 8.7$,
and $P<10^{-8}$ ($N$=234). At $0.1 < CR < 0.25$ cts/s,
$r_s = -0.27, T_S=-3.7$, and $P=0.00014$  ($N$=177), while at $CR>0.25$ cts/s, $r_s = -0.205, T_S = -2.7$, and $P = 0.0038$ ($N$=174), where $r_s$ is the Spearman rank order correlation coefficient,
$T_S$ is the Student T-test, $P$ is the probability that the distribution is just random, and $N$ is the number of data points.

\begin{table}
\scriptsize
        \centering
        \caption{
        X-ray spectral fit results of Mrk 335 (see Sect. 3.2).
        Columns are as follows: (1) Models: pl = powerlaw, bbdy = black body, zpcf = partial covering absorber, and zxipcf = ionized absorber model.
        (2) Absorption in units of 10$^{22}$ cm$^{-2}$; (3) power-law photon index; (4) $kT$ in units of keV of the black body, or ionization parameter log $\xi$ of the ionized absorber;
        (5) covering factor; and
        (6) goodness of fit, reduced $\chi^2$. If no errors are reported, the quantity was fixed. For all fits the column density at $z=0$ was fixed to the Galactic value of 3.96$\times 10^{20}$ cm$^{-2}$.}
        \label{tab:spec-fits}
        \begin{tabular}{lccccc}
        \noalign{\smallskip}
                \hline
                \noalign{\smallskip}
                model & $N_{\rm H}$ & $\Gamma_{\rm X}$ & $kT$ or log $\xi$ &
                $\eta$ &
                $\chi{^2}_{\rm red}$ \\
                (1) & (2) & (3) & (4) &  (5) & (6)   \\
                \noalign{\smallskip}
                \hline
                \noalign{\smallskip}
                \multicolumn{6}{c}{high state, 2020 June 23 -- July 1}  \\
                \noalign{\smallskip}
                \hline
                \noalign{\smallskip}
                pl & & 2.49\pl0.06 & &  & 3.2 \\
                zpcf * pl & 6.3$^{+3.3}_{-1.9}$ & 2.95$^{+0.14}_{-0.14}$ & &  0.86$^{+0.04}_{-0.05}$ &  1.4 \\
                bbdy+(zpcf * pl)$^1$ & 4.3$^{+1.7}_{-1.1}$ & 2.97\pl0.20 & 0.15 & 0.88$^{+0.04}_{-0.05}$ &  1.5\\
                zxipcf * pl & 16.1$^{+4.2}_{-4.9}$ & 2.3$^{+0.1}_{-0.1}$ & 2.1$^{+0.2}_{-0.1}$ & 0.87$^{+0.07}_{-0.06}$ & 1.1 \\
        \noalign{\smallskip}
                \hline
        \noalign{\smallskip}
                \multicolumn{6}{c}{deep low state, 2019 May$^2$}     \\
                \noalign{\smallskip}
                \hline
                \noalign{\smallskip}
                pl & & 2.34\pl0.20 & & & 1.6 \\
                zpcf * pl & 13.95$^{+8.62}_{-5.06}$ & 3.07\pl0.21 &  & 0.94$^{+0.03}_{-0.06}$ & 1.0  \\
                pl + bbdy & & 0.50$^{+0.51}_{-0.55}$ & 0.15\pl0.02 &  & 0.9  \\
                \noalign{\smallskip}
                \hline
        \end{tabular}

$^1$ bbdy parameters fixed at low-state value

$^2$ Fits using Cash statistics \citep{cash79}
\end{table}

\begin{figure}
\includegraphics[clip, trim=1.0cm 1.0cm 3.0cm 2.6cm, width=7.1cm]
{fig3}
    \caption{
    Selected Swift spectra of Mrk 335 in 2020 May--July.
    The 2020 June high state is plotted in black,
    an intermediate state from June 2 and 3 in red, and a low state from 2020
    (data from May 16, 23, and 27, and June 7) in green. All spectra are fit with
    a single power-law model with Galactic absorption, which displays strong
    fit residuals especially at low and high energies.}
    \label{fig:xray_spec}
\end{figure}

\subsection{Spectral energy distribution}
The SEDs during the high state in May-June 2020 in comparison with a low -state data set in May 2019 are shown in Fig. \ref{fig:SED}. The results
show that the optical-UV spectrum has become bluer. While the W2-U color
during the low state in 2019 was --0.03\pl0.06, it was --0.61\pl0.06 on 2020 May 31.
In X-rays,  Mrk 335 has not only become brighter, but the spectrum has
changed as well (see Sect. 3.2).

\section{Discussion}

The long-term Swift light curve of Mrk 335 reveals that its deepest
minimum-state in the UV was reached in 2018, when the X-ray emission was
also at a deep and long-lasting low state, with little variability left
in X-rays (Fig. \ref{fig:lc-Swift}).
Except for one deviating epoch in 2019, and with some short-time variability
superposed, the UV has been on a steady rise since its 2018 minimum,
reaching and exceeding its former Swift} peak value (2011) again in 2020 May-June
(Fig. \ref{fig:lc-Swift}){\footnote{at no time during the Swift monitoring
did the UV emission of Mrk 335 reach the highest brightness level recorded
with IUE in the 1980s (e.g., \citet{Dunn2006},
Fig. 1 of \citet{Grupe2008}, \citet{Longinotti2013}, Fig. 14 of \citet{Tripathi2020})}}.
Instead, the long-term X-ray light curve remained flat throughout 2018-2019,
and we only detected a significant rise since 2020 May, with multiple peaks.
Except that both reached a deep low-state in 2018 (\cite{Parker2019},
which was relatively short-lived in the UV, however; Fig. \ref{fig:lc-Swift}),
the variability observed with Swift in the UV and X-rays is not closely
correlated on the timescale of days to months \citep[see also][]{Grupe2012, Gallo2018, Tripathi2020}:
multiple flares in the UV are not seen in X-rays and vice versa, as
confirmed by our DCF analysis of the 2020 data. This behavior is markedly
different from that of other bright type 1
AGN  \citep{Shappee2014, Gardner2017, Buisson2017, McHardy2018, Edelson2019} and demands an explanation (see below).

\begin{figure}
\includegraphics[clip, trim=1.7cm 1.5cm 1.3cm 0.1cm, width=5.9cm, angle=270] {fig4a}
\includegraphics[clip, trim=1.7cm 1.5cm 1.3cm 0.1cm, width=5.9cm, angle=270] {fig4b}
    \caption{Hardness ratio vs. X-ray count rate for Swift XRT data of Mrk 335 since
2007 (black). Three regimes are apparent: a softer-when-fainter trend at lowest count rate, a softer-when-brigther
pattern at intermediate count rate, and a near-constant HR at highest count rate. The upper panel shows the full
count-rate range.
Green circles represent the 2020 outbursts. The lower panel zooms into the count-rate range up to 0.5 cts/s. The red
lines represent the fits to the three regimes (Sect. 3.3), and the blue squares mark the median in each bin.}
    \label{fig:HR}
\end{figure}
Several mechanisms are known to cause strong changes in the X-ray flux and spectral
states of AGN: The first mechanism is changes in ionization state or geometry of ionized or
neutral absorption that fully or partially covers the continuum source \citep[e.g.,][]{AbrassartCzerny2000, Turner2011, leighly2015}.
Varying absorption is frequently seen in intermediate-type Seyfert galaxies whose variability timescales can be as short as hours or days, but can also occur in
NLS1s \citep{Risaliti2011} when the absorber is part of a clumpy accretion-disk wind or
resides within the BLR.
The second mechanism is changes in the X-ray reflection of photons off the (inner) accretion disk
\citep{Ross2005}.
 This model has frequently been applied
to explain short-time spectral changes in type 1 Seyfert galaxies by changes in the
location and geometry of the corona.
Third, on longer timescales, changes in the accretion rate including disk instabilities
can cause high or low states in AGN
\citep[e.g.,][]{Ross2018, Czerny2019}. This last mechanism has been favored to explain
some of the extreme changing-look AGN that also switch their Seyfert types.

\begin{figure}
\includegraphics[clip, trim=1.9cm 1.5cm 2.0cm 0.8cm, width=6.0cm, angle=270]
{fig5}
    \caption{Observed Swift SED of Mrk 335 during 2020 June 2, June 27, and the deep low-state of 2019 May.
    Note the X-ray flux deficit around (0.7-1.2)\,keV on June 2. For comparison, the Swift SED during the highest X-ray state from 2009 January is shown (purple).
    }
    \label{fig:SED}
\end{figure}

In Mrk 335, the variability of the UV band is not correlated with the X-ray band on the
timescale of days to months \citep[with the exception of an epoch in 2014;][]{Gallo2018}.
This rejects models that predict a close link between UV and X-rays: For instance,
models where the observed X-rays are upscattered UV photons; or models where the UV
is reprocessed emission from an X-ray corona heating the disk, which have successfully been applied
to several well-observed AGN that show a close UV--X-ray
correlation \citep{Shappee2014, McHardy2018}; or models where UV and X-rays
are due to synchrotron radiation (unlikely for Mrk 335 even though it is a radio emitter).

Furthermore , we can also exclude models of distant {\em{dusty}} absorbers that diminish the
UV by extinction and simultaneously diminish the X-rays by absorption.
However, {\em{dust-free}} absorption, arising in the accretion-disk region for instance
in the form of a clumpy wind, is an efficient mechanism to cause X-ray spectral and flux
changes while leaving the broad-band UV unaffected. And our spectral fits of a
partial covering absorber consistently imply a decrease in covering factor
and column density as the source brightens in X-rays. Our fit parameters of
the ionization parameter, column density, and covering factor are consistent with predictions from models of clumpy disk winds \citep{Takeuchi2013}
that form via the Rayleigh-Taylor and radiation-hydrodynamic instability \citep{Shaviv2001} and
which predict variability timescales on the  order of weeks.

A partial covering scenario can also explain why Mrk 335 only shows a mild trend of softer-when-brighter (Fig. \ref{fig:HR}) because a partial coverer would preserve the spectral shape as it leaks different fractions of the same
intrinsic continuum when the covering factor changes. The softer-when-fainter pattern at the very lowest
count rates is likely due to a soft X-ray component in the form of multiple emission lines
\citep{Longinotti2019, Parker2019} and a (distant) reflection component \citep{Gallo2019} which become increasingly apparent as the primary continuum diminishes. While the X-ray covering factor obtained from the current and
previous (intermediate to low state) observations of Mrk 335 is high ($\sim$ 80-95\%), a partial covering scenario
for Mrk 335 is independently favored by previous UV observations that revealed variable absorption lines, including
Ly$\alpha$ at 20--30\%
covering fraction of the UV-emitting site \citep{Longinotti2019, Parker2019}.

There is other evidence that the X-rays we observe along our line of sight are not representative of the 4$\pi$ emission:
Optical reverberation mapping \citep{Grier2012} has shown that the Balmer line H$\beta$ and the
high-ionization line HeII$\lambda$4686 of Mrk 335 follow the amplitude of the optical continuum
variability{\footnote{the Balmer line emission is driven by the EUV beyond the Lyman limit, and the
ionization potential of HeII is 54 eV}} and not the higher amplitude seen in X-rays. While $R_{\rm max}$ = 1.53 in optical
continuum variability, and $R_{\rm max}$ = 1.55 in HeII (Grier et al. 2012; where $R_{\rm max}$ is the
ratio of maximum to minimum flux in each light curve), we measured during the same time interval $R_{\rm max}$ = 5.44 in X-rays.

Furthermore, no significant variability of the optical coronal lines in Mrk 335 was detected \citep{Grupe2008} despite
strong changes in observed X-rays. As those lines are driven by the EUV-to-soft-X-ray part of the SED, the
observations indicate that the emission lines saw a different (less variable) continuum.
An absorption scenario can also explain the lack of intrinsic variability of the photoionized soft X-ray emission lines
\citep{Parker2019}, and is consistent with the lack of IR variability of Mrk 335 \citep{Wright2010}.

The question remains why Mrk 335 is no optical Seyfert-type changer{\footnote{as observed in an increasing number of other AGN \citep{MacLeod2016}}}
despite its high amplitude of X-ray variability (a factor $\sim$50 between highest and lowest state in our
long-term Swift light curve). As described above, the emission-line regions seem to see a less variable EUV-X-ray
SED. Furthermore, because the broad Balmer lines are intrinsically bright, even a drop of a factor of 2 still preserves the
type 1 nature of Mrk 335.

As Mrk 335 is a very bright AGN at high state, it may turn out to be a Rosetta stone in deciphering the various
contributions of reflection, absorption, and intrinsic emission, which imprint their presence on AGN X-ray spectra.
Following Mrk 335 closely as it emerges from low state will therefore provide us with an important opportunity of
understanding the physical processes that drive long-term AGN state changes.

\begin{acknowledgements}
We would like to thank Brad Cenko for approving our Swift observing requests and the Swift team for carrying out
our observations, and our anonymous referee for very useful comments.
This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. ALL
acknowledges support from CONACyT grant CB-286316.
\end{acknowledgements}