Self-Supervised Single-Image Deconvolution with Siamese Neural Networks
Inverse problems in image reconstruction are fundamentally complicated by
unknown noise properties. Classical iterative deconvolution approaches amplify
noise and require careful parameter selection for an optimal trade-off between
sharpness and grain. Deep learning methods allow for flexible parametrization
of the noise and learning its properties directly from the data. Recently,
self-supervised blind-spot neural networks were successfully adopted for image
deconvolution by including a known point-spread function in the end-to-end
training. However, their practical application has been limited to 2D images in
the biomedical domain because it implies large kernels that are poorly
optimized. We tackle this problem with Fast Fourier Transform convolutions that
provide training speed-up in 3D microscopy deconvolution tasks. Further, we
propose to adopt a Siamese invariance loss for deconvolution and empirically
identify its optimal position in the neural network between blind-spot and full
image branches. The experimental results show that our improved framework
outperforms the previous state-of-the-art deconvolution methods with a known
point spread function.
Mikhail Papkov, Kaupo Palo, Leopold Parts
# Self-Supervised Single-Image Deconvolution with Siamese Neural Networks

###### Abstract

Inverse problems in image reconstruction are fundamentally complicated by unknown noise properties. Classical iterative deconvolution approaches amplify noise and require careful parameter selection for an optimal trade-off between sharpness and grain. Deep learning methods allow for flexible parametrization of the noise and learning its properties directly from the data. Recently, self-supervised blind-spot neural networks were successfully adopted for image deconvolution by including a known point-spread function in the end-to-end training. However, their practical application has been limited to 2D images in the biomedical domain because it implies large kernels that are poorly optimized. We tackle this problem with Fast Fourier Transform convolutions that provide training speed-up in 3D microscopy deconvolution tasks. Further, we propose to adopt a Siamese invariance loss for deconvolution and empirically identify its optimal position in the neural network between blind-spot and full image branches. The experimental results show that our improved framework outperforms the previous state-of-the-art deconvolution methods with a known point spread function.

Keywords:Deconvolution Microscopy Deep learning.

## 1 Introduction and related work

Quality enhancement is central for microscopy imaging [26, 29]. Its two most important steps are noise removal and blur reduction that improve performance in downstream tasks for both humans and algorithms.

Denoising, approaches to which were developed decades ago [5], was given a boost by recent advances in deep learning methods [1, 3, 10, 11, 14, 26, 15, 19, 18], which substantially outperform classical algorithms [5, 2]. A direct approach to train a denoising neural network is to provide it with pairs of low and high-quality images [26], but such data are rarely available in practice, which raises the need for alternative techniques. Lehtinen \(et\)\(al.\)[14] proposed to supervise the training with an independently acquired noisy copy of the input, while Batson and Royer [1] concurrently with Krull \(et\)\(al.\)[10] developed a theory of a \(J\)-invariant blind-spot denoising network that operates on a single image in a self-supervised manner.
The latter theory is based on the idea that a network cannot learn the noise that is conditionally pixel-wise independent given the signal, so learning to predict the values of masked pixels can improve the performance in restoration. Xie _et al_. [27] further questioned the necessity of blind spots for efficient denoising. They proposed Noise2Same with Siamese invariance loss between outputs of masked and unmasked inputs to prevent the network from learning an identity function. This advance pushed the performance of self-supervised denoising at a cost of doubled training time.

In addition to noise, equipment imperfection imposes blur on the images via a point spread function (PSF) [24]. The blur can be efficiently removed with classical iterative approaches such as Lucy-Richardson [21], but these algorithms tend to amplify the noise with each iteration and require careful regularization and stopping criteria [12]. As supervised methods are not applicable to deconvolution problem due to the ground truth inaccessibility, Lim _et al_. [16] proposed CycleGAN [30] with linear blur kernel, which performed well in both simulation and real microscopy deconvolution. However, the models of this class are notoriously hard to train [17] and tends to converge to a perceptually conceiving solution that does not necessarily reflect the true underlying structure [13]. Besides, it requires clean data examples, albeit unpaired. Deep self-supervised denoising systems were also adopted for deconvolution purposes. Kobayashi _et al_. [9] assumed an intermediate output of the network to be a deconvolved representation and proposed to train Noise2Self [1] as a pseudo-inverse model with a known PSF kernel. This method was only applied to 2D data.

Our contribution is three-fold. Firstly, we adopt the Siamese invariance loss [27] to the deconvolution task and identify the optimal neural network's outputs to apply it. Secondly, we train the first, to our knowledge, 3D self-supervised deconvolutional network without the adversarial component. Lastly, we propose to alleviate the computational costs of training a Siamese network by using the Fast Fourier Transform (FFT) for convolution with the PSF kernel.

## 2 Methods

We propose a generalized self-supervised deconvolution framework. We follow Kobayashi _et al_. [9] and define the framework as a composition of trainable model \(f(\cdot)\) followed by the fixed PSF convolution \(g(\cdot)\). The framework (Figure 1) allows for both one-pass blind-spot training[10, 1, 9] and two-pass training [27].

For the training, we use two inputs, unmasked \(\mathbf{x}\) and masked \(\mathbf{x}_{J^{c}}\) where the pixels at locations \(J\) are replaced with Gaussian noise. We optimize the trainable model using the composition loss \(\mathcal{L}(f)\) of six terms in Eq. (1).

**Blind-spot loss.**\(\mathcal{L}_{bsp}\) is a mean-squared error (MSE) between the network's unmasked input \(\mathbf{x}\) and output from the masked forward pass \(g(f(\mathbf{x}_{J^{c}}))_{J}\) measured at the locations of the masked pixels \(J\)[10, 1].

**Reconstruction loss.**\(\mathcal{L}_{rec}\) is an MSE between the network's unmasked input \(\mathbf{x}\) and output from the unmasked forward pass \(g(f(\mathbf{x}))\)[27].
**Reconvolved invariance loss.**\(\mathcal{L}_{inv}\) is an MSE between the network's output from the masked forward pass \(g(f(\mathbf{x}_{J^{c}}))_{J}\) and output from the unmasked forward pass \(g(f(\mathbf{x}))_{J}\).

**Deconvolved invariance loss.**\(\mathcal{L}_{inv\,(d)}\) is an MSE between the network's output from the masked forward pass \(f(\mathbf{x}_{J^{c}})_{J}\) and output from the unmasked forward pass \(f(\mathbf{x})_{J}\) before PSF convolution \(g(\cdot)\). Both Siamese invariance losses are computed between the network's outputs from altered (masked and unmasked) inputs at the location of the masked pixels \(J\). They prevent the network from learning an identity function from \(\mathcal{L}_{rec}\) minimization, which is especially important in the single-image deconvolution task.

**Boundary losses**. \(\mathcal{L}_{bound\,(d)}\) and \(\mathcal{L}_{bound}\) regularize the outputs to be within \([min,max]\) boundaries before or after the PSF convolution respectively [9]. These losses are measured after destandardization. We use \(min=0\) and \(max=1\).

In this work, we compare three practical cases. The first one, Eq. (2), considers only the blind-spot loss \(\mathcal{L}_{bsp}\) and establishes a self-supervised \(J\)-invariant baseline.

Figure 1: The self-supervised deconvolution framework. The trainable model \(f(\cdot)\) takes noisy input \(\mathbf{x}\) in unmasked forward pass (solid) and masked noisy input \(\mathbf{x}_{J^{c}}\) in masked forward pass (dashed). The model produces two intermediate outputs (deconvolved representations). These outputs are convolved with a fixed PSF convolution \(g(\cdot)\) to obtain final outputs (reconvolved representations). The deconvolved and reconvolved invariance losses are computed between the respective outputs before and after the PSF. The reconstruction loss is computed between the noisy input and the unmasked forward pass output. The blind-spot loss is computed between the noisy input and masked forward pass output, it does not require an unmasked forward pass.
The second and the third ones are not \(J\)-invariant, both of them include reconstruction loss and invariance loss calculated before (Eq. (4)) or after (Eq. (3)) the PSF convolution. By default, we calculate boundary regularization loss and invariance loss from the same output. We set \(\lambda_{inv}\) and \(\lambda_{inv\,(d)}\) to 2 [27], \(\lambda_{bound}\) and \(\lambda_{bound\,(d)}\) to 0.1 [9] unless otherwise specified.

\[\begin{split}\mathcal{L}(f)=&\lambda_{bsp}\, \mathbb{E}_{J}\mathbb{E}_{x}\|g(f(\mathbf{x}_{J^{c}}))_{J}-\mathbf{x}_{J}\|^{2 }+\\ &\lambda_{rec}\,\mathbb{E}_{x}\|g(f(\mathbf{x}))-\mathbf{x}\|^{2} /m+\\ &\lambda_{inv}\,\mathbb{E}_{J}\left[\mathbb{E}_{x}\|g(f(\mathbf{x }))_{J}-g(f(\mathbf{x}_{J^{c}}))_{J}\|^{2}/|J|\right]^{1/2}+\\ &\lambda_{inv\,(d)}\,\mathbb{E}_{J}\left[\mathbb{E}_{x}\|f( \mathbf{x})_{J}-f(\mathbf{x}_{J^{c}})_{J}\|^{2}/|J|\right]^{1/2}+\\ &\lambda_{bound}\,\mathbb{E}_{x}\left(|min-g(f(\mathbf{x}))|+|g( f(\mathbf{x}))-max|\right)+\\ &\lambda_{bound\,(d)}\,\mathbb{E}_{x}\left(|min-f(\mathbf{x})|+|f( \mathbf{x})-max|\right)\end{split} \tag{1}\]

\[\mathcal{L}(f)_{\text{Noise2Self}} =\mathcal{L}_{bsp} \tag{2}\] \[\mathcal{L}(f)_{\text{Noise2Same}} =\mathcal{L}_{rec}+\lambda_{inv}\,\mathcal{L}_{inv}+\lambda_{bound }\,\mathcal{L}_{bound}\] (3) \[\mathcal{L}(f)_{\text{Noise2Same\,(d)}} =\mathcal{L}_{rec}+\lambda_{inv\,(d)}\,\mathcal{L}_{inv\,(d)}+ \lambda_{bound\,(d)}\,\mathcal{L}_{bound\,(d)} \tag{4}\]

#### 3.2.2 Neural network architecture

For 2D data, we use the previously proposed [10, 27] variation of a U-Net [22] architecture without modifications. The fully-convolutional network of depth 3 consists of an encoder and a decoder with concatenating skip connections at corresponding levels. The first convolutional layer outputs 96 features, and this number doubles with every downsampling step. For 3D data, we modify the network to have 48 output features from the first convolution. We also replace the concatenation operation with an addition in skip connections. Such modification allows for a reduction in training time for 3D convolutional networks without substantial performance sacrifices. We implement the network using PyTorch [20].

#### 3.2.3 Training

During training, we sample a patch (\(128\times 128\) pixels for 2D images, \(64\times 64\times 64\) for 3D), randomly rotate and flip it. We use batch sizes of 16 for 2D images and 4 for 3D images. Each image is standardized individually to zero mean and unit variance. For the masked forward pass, we randomly replace 0.5% of pixels from each training image with Gaussian noise (\(\sigma=0.2\)) [27]. We use Adam [8] optimizer and multiply its learning rate from \(4\cdot 10^{-4}\) by 0.5 every 500 steps out of 3k total for 2D data [9] and every 2k steps out of 15k total for 3D data. We train and evaluate all models on a single NVIDIA V100 16GB GPU.

#### 3.2.4 Inference

For inference, we do only unmasked forward pass and discard the PSF convolution during inference [9]. We use the last model checkpoint since we do not possess any reliable validation metric. For large images that do not fit in
GPU memory, we predict by overlapping patches and stitch the output together using pyramid weights [7] which allows for avoiding edge artifacts. In 3D data, we predict patches of \(128\times 128\times 128\) with overlaps of 32 pixels.

#### 2.0.2 Fast Fourier Transform Convolution

We use FFT for the convolution with a fixed PSF during training. Convolution operation \(x*k\) is equivalent to element-wise matrix multiplication in Fourier space \(F^{-1}(F(x)\odot F(k))\). If image \(x\) is \(M\times M\) pixels and kernel \(k\) is \(N\times N\) and both of them have dimensionality of \(d\), the complexity of an ordinary convolution is \(O(M^{d}N^{d})\), while the complexity of a Fourier convolution is \(O(M^{d}d\log M)\) and does not depend on \(N\). With large enough \(N\), convolution in Fourier space becomes computationally cheaper, and large kernels (even of the size of the image [24]) are common for PSF. It is recommended to use FFT with 2D kernels for \(N>25\) and with 3D kernels for \(N>9\)[6].

## 3 Experiments

We validate our approach on 2D and 3D data in a single image deconvolution task. In all the experiments we train one model per image to obtain a lower bound for the extreme self-supervised case. Every image is normalized within \([0,1]\). We generate blurry images by convolving them with a realistic Richards & Wolf PSF of a 0.8NA 16x microscope objective with \(17\times 17\) or \(17\times 17\times 17\) 0.406-micron pixels. Then we add a mixture of Poisson (\(\alpha=0.001\)), Gaussian (\(\sigma=0.1\)), and salt-and-pepper (only for 2D images, \(p=0.01\)) noise to these images and quantize them with 10 bit precision [9].

We compare the results against classical Lucy-Richardson (LR) deconvolution algorithm [21] and deep learning Self-Supervised Inversion (SSI) [9] algorithm. For LR evaluate the results after 2, 5, 10, and 20 iterations to observe the trade-off between deconvolution quality and noise amplification and report the best ones. We do not include in comparison other classical methods such as Conjugate Gradient optimization [4] and Chambole-Pock primal-dual inversion [4], because they were shown inferior to both LR and SSI [9]. We rerun the baseline experiments ourselves wherever possible.

We evaluate the model performance against clean images using root mean squared error (RMSE), peak signal-to-noise ratio (PSNR) [28], and structural similarity index(SSIM) [25]. For PSNR and SSIM, we explicitly set the data range to the true values \([0,1]\). Since it was not done in prior work [9], some values in Table 2 are missing to avoid confusion. Additionally, for 2D images, we report mutual information (MI) [23] and spectral mutual information (SMI) [9]. For all metrics except RMSE, higher is better.

### 2D dataset

First, we tested the deconvolution framework performance for the 2D case on a benchmark dataset of 22 single-channeled images [9] with size varying from
\(512\times 512\) to \(2592\times 1728\) pixels. We compared three different variants of loss function, Eqs. (2)- (4) (Table 1). Loss (4) performed best with \(\lambda_{inv\,(d)}=2\) and \(\lambda_{bound\,(d)}=0.1\), achieving \(\text{PSNR}=22.79\), \(\text{SSIM}=0.46\), and \(\text{RMSE}=0.078\).

We compared our model against the baseline LR and SSI methods (Table 2, Figure 2). It performed best by PSNR, SMI, and RMSE, and showed similar results by SSIM (0.01 less than LR) and MI (0.01 less than SSI).

### 3D dataset

We evaluated the framework performance for the 3D case on a single synthetic image of microtubules [24] of size \(128\times 256\times 512\). From the three variants of loss function, Eqs. (2)- (4), loss (3) performed best with \(\lambda_{inv}=2\) and \(\lambda_{bound}=0\) showing the results \(\text{PSNR}=24.08\), \(\text{SSIM}=0.39\), and \(\text{RMSE}=0.063\) (Table 3).

We also compared the best-performing model against the LR and SSI methods (Table 4, Figure 2). Our model surpassed baselines by a substantial margin (+1.5 PSNR against SSI and +1.7 against LR).

## 4 Discussion

Our Siamese neural network was superior in single-image deconvolution against both classical and deep learning baselines. The architecture is simple and did not require additional tricks such as masking schedule or adding noise to gradients [9] for training. Training converged to similar performance for all of the several tested random weight initializations.

Loss function performance was inconsistent between the 2D and 3D cases. Invariance loss applied to deconvolved representation proved best for 2D data but

Figure 2: Example deconvolution results on _Drosophila_ 2D image (top), microtubules 3D image, frontal max projection (bottom), magnified. For Lucy Richardson, the result is shown after the best iteration. Our algorithm uses loss from Eq. (4) with \(\lambda_{inv\,(d)}=2\), \(\lambda_{bound\,(d)}=0.1\). PSNR metric is reported for the presented images.
failed to give an advantage in 3D: despite the high SSIM, images appear noisy. We hypothesize that the optimal solution depends on the data; _e.g._ in microtubules, the 3D dataset signal is very sparse. Boundary loss did not drastically affect the training of Siamese networks.

Despite performance superiority, Siamese networks are twice more expensive in computation because they require two forward passes through the neural network. Additionally, Noise2Same U-Net architecture has 10x more trainable parameters than SSI (5.75M against 0.55M in 2D). This problem is exacerbated by the necessity of convolutions with large PSF kernels which are not optimized in modern GPUs. We propose to alleviate this problem by using FFT for convolution with PSF. It leads to 2.3x speed improvement for 3D deconvolution with PSF of size \(17\times 17\times 17\), and this advantage will grow for larger kernels. For example, while training with a PSF of \(31\times 31\times 31\) we observed a 14x speedup.

\begin{table}
\begin{tabular}{l c c c c c c c c} \hline \hline  & \multicolumn{3}{c}{\(\lambda_{bsp}\)\(\lambda_{rec}\)\(\lambda_{inv\,(d)}\)\(\lambda_{inv}\)\(\lambda_{bound\,(d)}\)\(\lambda_{bound}\) PSNR \(\uparrow\) SSIM \(\uparrow\) RMSE \(\downarrow\)} \\ \hline input & & & & & & 17.77 & 0.18 & 0.130 \\ \hline \(\mathcal{L}(f)_{\text{Noise2Self}}\) & 1 & 0 & 0 & 0 & 0 & 0 & 22.49 & 0.27 & 0.079 \\  & 1 & 0 & 0 & 0 & 0 & 0.1 & 22.51 & 0.43 & 0.079 \\ \hline \(\mathcal{L}(f)_{\text{Noise2Same}}\) & 0 & 1 & 0 & 2 & 0 & 0 & 21.49 & 0.23 & 0.089 \\  & 0 & 1 & 0 & 2 & 0 & 0.1 & 21.81 & 0.41 & 0.085 \\ \hline \(\mathcal{L}(f)_{\text{Noise2Same (d)}}\) & 0 & 1 & 2 & 0 & 0 & 0 & **22.79** & 0.46 & **0.078** \\  & 0 & 1 & 2 & 0 & 0.1 & 0 & 22.67 & **0.46** & 0.078 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Denoising results for the 2D dataset with different coefficients for loss (1) components: \(\lambda_{bsp}\) — coefficient for the blind-spot loss, \(\lambda_{rec}\) — coefficient for the reconstruction loss, \(\lambda_{inv\,(d)}\) — coefficient for deconvolved invariance loss, \(\lambda_{inv}\) — coefficient for reconvolved invariance loss, \(\lambda_{bound}\) and \(\lambda_{bound\,(d)}\) — coefficients for boundary regularization loss for reconvolved and deconvolved images. For all metrics except RMSE, higher is better. The best values are highlighted in bold.

\begin{table}
\begin{tabular}{l c c c c c c c} \hline \hline  & \multicolumn{3}{c}{PSNR \(\uparrow\) SSIM \(\uparrow\) MI \(\uparrow\) SMI \(\uparrow\) RMSE \(\downarrow\) Train \(t\) (s) Inference \(t\) (ms)} \\ \hline input & 17.8 & 0.18 & 0.07 & 0.18 & 0.131 & - & - \\ \hline LR \(n=2\) & 22.0 & **0.47** & 0.13 & 0.25 & 0.086 & - & 46 \\ LR \(n=5\) & 22.2 & 0.44 & 0.12 & 0.25 & 0.082 & - & 48 \\ SSI (reported [9]) & 22.5 & - & **0.14** & 0.27 & - & 238 & 26 \\ SSI (reproduced) & 22.0 & 0.46 & 0.14 & 0.26 & - & 442 & 26 \\ ours & **22.8** & 0.46 & 0.13 & **0.30** & **0.078** & 856 & 53 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Our best-performing method for 2D data compared to Lucy Richardson (LR) [21] and Self-Supervised Inversion (SSI) [9] baselines by denoising metrics (see Section 3 for details), training and inference time (measured for _Drosophila_ image \(1352\times 532\) pixels). For all metrics except RMSE, higher is better. The best values are highlighted in bold. Our algorithm uses loss (4) with \(\lambda_{inv\,(d)}=2\), \(\lambda_{bound\,(d)}=0.1\)
## 5 Conclusion

Blind-spot networks [10, 1] were seminal in self-supervised denoising and performed comparably to supervised methods [14, 26]. Their success was later translated to image deconvolution [9], a more complicated inverse problem. However, it was shown that \(J\)-invariance leads to suboptimal performance in denoising because masked pixels contain useful bits of information [27]. In this work, we presented a novel unified image deconvolution framework, which generalized the accumulated prior advances, and set a new standard for image quality enhancement performance. We investigated the contributions of various self-supervised loss components and empirically identified the optimal usage scenarios for 2D and 3D data. We also proposed using Fast Fourier Transform for the training of deconvolution neural networks, which drastically speeds up the computation, especially for large kernels. This advantage allows one to use our method for non-confocal microscopy with extremely big point spread functions.

\begin{table}
\begin{tabular}{l c c c c c c} \hline \hline  & \multicolumn{2}{c}{PSNR \(\uparrow\) SSIM \(\uparrow\) RMSE \(\downarrow\) Train \(t\) (s) Inference \(t\) (ms)} \\ \hline input & 21.2 & 0.64 & 0.087 & - & - \\ \hline LR \(n=2\) & 22.4 & 0.32 & 0.076 & - & 69 \\ SSI & 22.6 & 0.39 & 0.074 & 2814 (6580) & 578 \\ ours & **24.1** & **0.39** & **0.062** & 8160 (18900) & 734 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Our best performing method for 3D data compared to Lucy Richardson (LR) [21] and Self-Supervised Inversion (SSI) [9] baselines by denoising metrics (see Section 3 for details), training and inference time. For all metrics except RMSE, higher is better. The best values are highlighted in bold. Our algorithm uses loss (3) with \(\lambda_{inv}=2\), \(\lambda_{bound}=0\). Training time is reported with and without FFT (in brackets)

\begin{table}
\begin{tabular}{l c c c c c c c c} \hline \hline  & \multicolumn{2}{c}{\(\lambda_{bsp}\)\(\lambda_{rec}\)\(\lambda_{inv\,(d)}\)\(\lambda_{inv}\)\(\lambda_{bound\,(d)}\)\(\lambda_{bound}\) PSNR \(\uparrow\) SSIM \(\uparrow\) RMSE \(\downarrow\)} \\ \hline input & & & & & & 21.21 & 0.64 & 0.087 \\ \hline \(\mathcal{L}(f)_{\text{Noise2Self}}\) & 1 & 0 & 0 & 0 & 0 & 0 & 23.71 & 0.39 & 0.065 \\  & 1 & 0 & 0 & 0 & 0 & 0.1 & 23.71 & 0.39 & 0.065 \\ \hline \(\mathcal{L}(f)_{\text{Noise2Same}}\) & 0 & 1 & 0 & 2 & 0 & 0 & **24.08** & 0.39 & **0.063** \\  & 0 & 1 & 0 & 2 & 0 & 0.1 & 24.02 & 0.39 & 0.063 \\ \hline \(\mathcal{L}(f)_{\text{Noise2Same (d)}}\) & 0 & 1 & 2 & 0 & 0 & 0 & 21.05 & **0.82** & 0.089 \\  & 0 & 1 & 2 & 0 & 0.1 & 0 & 21.05 & **0.82** & 0.089 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Denoising results for the 3D dataset with different coefficients for loss (1) components: \(\lambda_{bsp}\) — coefficient for the blind-spot loss, \(\lambda_{rec}\) — coefficient for the reconstruction loss, \(\lambda_{inv\,(d)}\) — coefficient for deconvolved invariance loss, \(\lambda_{inv}\) — coefficient for reconvolved invariance loss, \(\lambda_{bound}\) and \(\lambda_{bound\,(d)}\) — coefficients for boundary regularization loss for reconvolved and deconvolved images. For all metrics except RMSE, higher is better. The best values are highlighted in bold.
## Acknowledgements

This work was funded by Revvity Inc. (VLTAT19682) and Wellcome Trust (206194). We thank High Performance Computing Center of the Institute of Computer Science at the University of Tartu for the provided computing power.
