Multi-scale Conditional Generative Modeling for Microscopic Image Restoration

## Abstract

Abstract The advance of diffusion-based generative models in recent years has revolutionized state-of-the-art (SOTA) techniques in a wide variety of image analysis and synthesis tasks, whereas their adaptation on image restoration, particularly within computational microscopy remains theoretically and empirically underexplored. In this research, we introduce a multi-scale generative model that enhances conditional image restoration through a novel exploitation of the Brownian Bridge process within wavelet domain. By initiating the Brownian Bridge diffusion process specifically at the lowest-frequency subband and applying generative adversarial networks at subsequent multi-scale high-frequency subbands in the wavelet domain, our method provides significant acceleration during training and sampling while sustaining a high image generation quality and diversity on par with SOTA diffusion models. Experimental results on various computational microscopy and imaging tasks confirm our method’s robust performance and its considerable reduction in its sampling steps and time. This pioneering technique offers an efficient image restoration framework that harmonizes efficiency with quality, signifying a major stride in incorporating cutting-edge generative models into computational microscopy workflows.

## 1 Introduction

Within the last decade, the landscape of image synthesis has been radically transformed by the advent of generative models (GMs) [, , ]. Among their broad success in various image synthesis applications, image restoration, including super-resolution, shadow removal, inpainting, etc, have caught much attention due to their importance in various practical scenarios. Image restoration aims to recover high-quality target image from low-quality images measured by an imaging system with assorted degradation effects, e.g., downsampling, aberration and noise. Numerous tasks in photography, sensing and microscopy can be formulated as image restoration problems, and therefore the importance of image restoration algorithms is self-evident in practical scenarios [, , , ].

Due to the ill-posedness of most image restoration problems, the application of generative learning becomes crucial for achieving high-quality image reconstruction. The wide applications of deep learning (DL)-based generative models in image restoration began with the success of generative adversarial networks (GANs) []. The emergence of conditional GANs achieved unprecedented success and outperformed classical algorithms in various applications of computational imaging and microscopy [, , , , , , , ]. Although Generative Adversarial Networks (GANs) have achieved remarkable success in image restoration, they are also known to be prone to training instability and mode collapse. These issues significantly restrict the diversity of outputs produced by GAN models [, ]. Recently, diffusion models (DMs) were intensively studied and outperformed GANs in various image generation tasks [, ]. Derived from the stochastic diffusion process, DMs employ neural networks to approximate the reverse diffusion process and sample target images from white noise with high quality and good mode coverage. Moreover, DMs have also been applied to image restoration, including super-resolution, dehazing, colorization, etc [, , , , , ]. Most of the existing methods regard the low-quality image as one additional condition in the reverse process and utilize it as one of the network arguments for inference. However, these image restoration models lack a clear modelling of the conditional image in the forward process and a theoretical guarantee that the terminal state of the diffusion is closely related to and dependent on the condition.

Figure: Figure 1: Inherent information loss in autoencoder backbones. The same high-quality image passes the three generative models without going through the diffusion process. Public pre-trained checkpoints are used for LDM [] and Refusion []. LDM and Refusion suffer from reconstruction loss due to lossy autoencoder backbones while the wavelet and inverse wavelet transform pairs are lossless.
Refer to caption: x1.png

Besides, compared to GANs, DMs are further limited by time-consuming sampling and iterative refinement. Many efforts have been made to improve the sampling speed while ensuring the quality of image generation, among which the Latent Diffusion Model (LDM) [] stands out as a notable example. LDM projects the diffusion process into a low-dimensional latent space of a pre-trained autoencoder, and hence greatly reduces the computational resource required for high-resolution image generation and restoration. The success of LDM inspired a few studies leveraging the latent representations of pretrained networks to reconstruct high-quality images from low-quality images [, ]. However, the performance of diffusion models in these approaches is highly dependent on the quality of latent representations of their pretrained backbones. There exist inevitable trade-offs between the compactness of the latent space and the reconstruction fidelity of the autoencoder, and between the normality of its latent distribution and its sample diversity [, , ]. Some study, e.g., Refusion, has explored to improve the fidelity of autoencoder by adding hidden connections between the encoder and decoder [], whereas the inherent information loss of learned compression cannot be completely eliminated. Though these methods reported improved perceptual scores between their outputs and reference images, in image restoration tasks where pixel-wise structural correspondence is emphasized, the reconstruction quality is constrained by the fidelity of autoencoder backbones, producing noticeable distortions and hallucinations as shown in Figure 1.

On the other hand, microscopy image restoration problems are usually composed of complicated degradation models, sparse image distributions and strict constraints on structural correspondence of the outputs to the conditional images [, , ]. To preserve the signal sparsity of common biological samples and local feature correspondence to the conditional image, existing methods leveraged CNN-based models and trained on carefully registered datasets [, , ]. Moreover, to enhance training and inference stability, conditional GAN were widely utilized to directly predict the corresponding ground truth image in a deterministic manner []. However, the adaptability and effectiveness of advanced generative models like DMs to microscopy image restoration tasks are under-explored.

To address the aforementioned limitations, here we introduce a novel multi-scale conditional generative model (MSCGM) for image restoration based on multi-scale wavelet transform and Brownian bridge stochastic process. For one thing, multi-scale wavelet transform effectively and losslessly compresses the spatial dimensions of conditional images, eliminating the lossy encoding process of current autoencoder-based DMs. Notably, our method does not involve the pre-training of autoencoders in competitive methods on a sufficiently large and diverse dataset sampled from the image domain. For another, Brownian bridge stochastic process incorporates the modelling of low-quality conditional image into both the forward and reverse diffusion process and better utilizes the information of the conditional images. Besides, we theoretically analyze the distributions of low- and high-frequency wavelet subbands and apply Brownian bridge diffusion process (BBDP) and GAN to the multi-scale generation of low- and high-frequency subbands, respectively. In sum, the contributions of this work are three-fold:

1. We are the first to factorize the conditional image generation within multi-scale wavelet domains, establishing a theoretical groundwork for a multi-scale conditional generative modeling;

2. Capitalizing on the unique distribution characteristics of wavelet subbands, we propose the innovative MSCGM, which seamlessly integrates BBDP and GAN;

3. We evaluate the MSCGM on various image restoration tasks, demonstrating its superior performance in both sampling speed and image quality than competitive methods.

## 2 Related works

### 2.1 Microscopy image restoration

As an important branch of computational imaging, computational microscopy springs up in recent years and aims to restore high-quality, multi-dimensional images from low-quality, low-dimensional measurements, usually with under-resourced equipment. Since the first work on microscopy image super-resolution reported in 2017 [], DL has enabled a wide spectrum of novel applications that were impossible with conventional optical technologies, e.g., microscopy image super-resolution surpassing the physical resolution limit of microscopic imaging systems [], volumetric imaging reconstructing 3D sample volumes from sparse 2D measurements [], and virtual image labelling to match the contrast conventionally provided by chemical or biological markers []. Compared to general image restoration in computational imaging, microscopy image restoration mainly differs in two aspects: (1) The degeneration process, including the transfer function, noise and aberration of the imaging system is generally complex, unknown and hard to measure precisely; and such degeneration process could vary significantly in real-world scenarios due to the variations of subjects, hardware and imaging protocols. (2) Strict pixel-wise correspondence between output and ground truth images and consistency with physical laws are generally emphasized [].

### 2.2 Generative models

GANs are well-known for generating high-qaulity, photorealistic samples rapidly [, ]. Through training a discriminator that tells ground truth images apart from fake images generated by the generator network, GANs outperformed traditional CNNs trained with hand-crafted, pixel-based structural losses such as $L_{1}$ and $L_{2}$ by providing a high-level, learnable perceptual objective. Conditional GANs such as Pix2Pix [], pix2pix HD [] and starGAN [] have been successfully applied in a wide spectrum of image-to-image translation and image restoration tasks, including image colorization [], style transfer [], image deblurring [], etc. Unsupervised image-to-image translation has also been extensively explored, such as cycleGAN [], UNIT [], DualGAN [], etc. In the fields of biomedical and microscopy imaging, researchers have also explored the applications of GANs, e.g., reconstructing low-dose CT and MRI images [, ], denoising microscopy images [], super-resolving diffraction-limited microscopy images [], among others.

Alternatively, transformer and related architectures have recently emerged and shown superior performance over convolutional neural networks (CNNs)-based GANs. Swin transformer [] and SwinIR [] have established a strong transformer-based baseline with competitive performance to CNNs. TransGAN substituted CNNs in the common GAN framework with transformers and improved the overall performance []. More recently, diffusion models (DMs) have been introduced and proved to be the state-of-the-art generative model in various image generation benchmarks [, ]. Dhariwal et al. reported diffusion models beat GANs on various image synthesis tasks []. Nevertheless, the democratization of diffusion models was limited by its huge demand for computational resources in both training and sampling. Recent advancements in fast sampling diffusion models [, , , , , ] have accelerated the sampling process in image generation tasks. DDGAN [] effectively combines the strengths of GANs with diffusion principles for quicker outputs, while LDM utilizes a compressed latent space to speed up generation. DDIM [] optimizes denoising steps for efficiency, FastDPM [] focuses on algorithmic improvements for rapid sampling, and CDM [] employs a staged approach to simplify the diffusion process.

## 3 Methods

### 3.1 Preliminaries

The outstanding success and wide applications of score-based diffusion models have been witnessed in the past years. Generally, for a Gaussian process $\{{\bm{x}}_{t},t=1,\cdots,T\}$ defined as:

$$ $q({\bm{x}}_{t}|{\bm{x}}_{t-1})=\mathscr{N}({\bm{x}}_{t};\sqrt{1-\beta_{t}}{\bm {x}}_{t-1},\beta_{t}{\bm{I}}),t=1,\cdots,T,$ (1) $$

the denoising DM attempts to solve the reverse process by parameterizing the conditional reverse distribution as:

$$ $p_{\theta}({\bm{x}}_{t-1}|{\bm{x}}_{t})=\mathscr{N}\Big{(}{\bm{x}}_{t};\frac{1 }{\sqrt{\alpha_{t}}}\big{(}{\bm{x}}_{t}-\frac{1-\alpha_{t}}{\sqrt{1-\bar{ \alpha}_{t}}}\bm{\epsilon}_{\theta}({\bm{x}}_{t},t)\big{)},\sigma_{t}{\bm{I}} \Big{)}.$ (2) $$

Here $\alpha_{t},\bar{\alpha}_{t},\beta_{t},\sigma_{t}$ are constants, and $\bm{\epsilon}_{\theta}$ is the network estimating the mean value of the reverse process. Despite the high image quality achieved by the denoising process, the total sampling steps $T$ can be very large and make the sampling process time-consuming.

###### Theorem 1 .

For a given error $\varepsilon$ between the generated distribution $p_{\theta}$ and the true distribution $p$, the sampling steps needed can be expressed by []:

$$ $T=O(\varepsilon^{-2}\kappa^{3}),$ (3) $$

where $\kappa$ is the condition number of the covariance matrix of $p$.

The proof of Theorem 1 is detailed in Appendix B.

Therefore, for highly non-Gaussian distributed images, e.g., microscopy images (please refer to Appendix C, where we theoretically and statistically demonstrate the high degree of non-Gaussianity of microscopy datasets), which tends to have high sparsity, resolution and contrast, standard diffusion models may not be practical due to their slow sampling speed and large sampling steps required for such distributions.

To overcome this limitation, we turned to multi-scale wavelet transform (as detailed in Appendix D), which offers an excellent latent space (i.e., wavelet domain) for generative modeling. The wavelet domain not only facilitates lossless compression but also provides low-frequency coefficients with near-Gaussian distributions. Consequently, we transform the diffusion process into the wavelet domain, thereby introducing a novel multi-scale wavelet-based generative model. It can be shown that the diffusion in wavelet domain is a dual problem to the original diffusion in the spatial domain. For a detailed explanation of this duality, please see Appendix E. Due to the distinct characteristics of wavelet coefficients in low- and high-frequency subbands (as detailed in Appendix C), we adopt different generative modeling approaches for the low- and high-frequency coefficients, marking a key innovation in our work. Specifically, while the low-frequency wavelet coefficients exhibit a Gaussian tendency, the high-frequency coefficients are sparse and non-Gaussian. Therefore, for the low-frequency coefficients, we employed the Brownian Bridge Diffusion Process (BBDP) [], and for the high-frequency coefficients, we utilized a Generative Adversarial Network (GAN) based generative method.

###### Theorem 1 .

### 3.2 Brownian bridge diffusion process

We leverage BBDP to better model the conditional diffusion process and apply it to image restoration. Image restoration tasks focus on the generation of the target image ${\bm{x}}_{0}\in\mathbb{R}^{H\times W\times C}$ from a conditional image ${\bm{y}}\in\mathbb{R}^{H\times W\times C}$. Most existing DMs tackle the conditional image ${\bm{y}}$ as an additional input argument to $\bm{\epsilon}_{\theta}$, without integrating the conditional probability distribution on ${\bm{y}}$ into the diffusion theory. Different from the standard diffusion process, we adapt the Brownian bridge process and derive a conditional diffusion process for image-to-image translation, termed as BBDP.

###### Definition 2 (Forward Brownian bridge process) .

The forward Brownian bridge with initial state ${\bm{x}}_{0}$ and terminal state ${\bm{y}}$ is defined as

$$ $q({\bm{x}}_{t}|{\bm{x}}_{0},{\bm{x}}_{T}={\bm{y}})=\mathscr{N}\Big{(}(1-\frac{ t}{T}){\bm{x}}_{0}+\frac{t}{T}{\bm{y}},\frac{t(T-t)}{T^{2}}{\bm{I}}\Big{)},$ (4) $$

By denoting $m_{t}=\frac{t}{T}$ and $\delta_{t}=\frac{t(T-t)}{T^{2}}$, we can reparameterize the distribution of ${\bm{x}}_{t}$ as:

$$ ${\bm{x}}_{t}={\bm{x}}_{0}+m_{t}({\bm{y}}-{\bm{x}}_{0})+\sqrt{\delta_{t}}\bm{ \epsilon}_{t},\\ \bm{\epsilon}_{t}\sim N(\bm{0},{\bm{I}})$ (5) $$

Following existing works on diffusion models, a network is trained to estimate ${\bm{x}}_{0}$ from ${\bm{x}}_{t},{\bm{y}}$; in other words, a network $\bm{\epsilon}_{\theta}$ is trained to estimate $m_{t}({\bm{y}}-{\bm{x}}_{0})+\sqrt{\delta_{t}}\bm{\epsilon}_{t}$.

The loss function of $\bm{\epsilon}_{\theta}$ is defined as:

$$ $L=\Sigma_{t}\gamma_{t}\mathbb{E}_{({\bm{x}}_{0},{\bm{y}}),\bm{\epsilon}_{t}} \lVert m_{t}({\bm{y}}-{\bm{x}}_{0})+\sqrt{\delta_{t}}\bm{\epsilon}_{t}-\bm{ \epsilon}_{\theta}({\bm{x}}_{t},{\bm{y}},t)\rVert_{2}^{2},$ (6) $$

where $\gamma_{t}$ is the weight for each $t$.

###### Theorem 3 .

The reverse process can be shown to be a Gaussian process with mean $\bm{\mu}^{\prime}_{t}({\bm{x}}_{t},{\bm{x}}_{0},{\bm{y}})$ and variance $\delta_{t}^{\prime}{\bm{I}}$

$$ $\displaystyle\begin{aligned} p({\bm{x}}_{t-1}|{\bm{x}}_{t},{\bm{x}}_{0},{\bm{y }})=\mathscr{N}({\bm{x}}_{t-1};\bm{\mu}^{\prime}_{t}({\bm{x}}_{t},{\bm{x}}_{0} ,{\bm{y}}),\delta_{t}^{\prime}{\bm{I}}),\end{aligned}$ (7) $\displaystyle\begin{aligned} \bm{\mu}^{\prime}_{t}({\bm{x}}_{t},{\bm{x}}_{0},{ \bm{y}})=c_{xt}{\bm{x}}_{t}+c_{yt}{\bm{y}}-c_{\epsilon t}\bm{\epsilon}_{\theta }({\bm{x}}_{t},{\bm{y}},t)\end{aligned}$ (8) $\displaystyle\delta_{t}^{\prime}=\frac{\delta_{t|t-1}\delta_{t-1}}{\delta_{t}}.$ (9) $$

Please refer to Appendix F regarding the expression of $c_{xt},c_{yt},c_{\epsilon t}$ and a detailed proof. Equation 8 indicates the posterior sampling process of BBDP and the training process of BBDP is detailed in Appendix Algorithm 1.

###### Definition 2 (Forward Brownian bridge process) .

###### Theorem 3 .

### 3.3 Multi-scale conditional generative model

Wavelet transform, with its theoretical details outlined in Appendix D, is characterized by an orthogonal transform matrix ${\bm{A}}\in\mathbb{R}^{N^{2}\times N^{2}}$. The wavelet transform decomposes an image ${\bm{x}}\in\mathbb{R}^{N^{2}}$ to one low-frequency (LL) subband ${\bm{x}}^{1}_{L}\in\mathbb{R}^{\frac{N^{2}}{4}}$ and remaining high-frequency subbands ${\bm{x}}^{1}_{H}\in\mathbb{R}^{\frac{3N^{2}}{4}}$.

###### Definition 4 (Multi-scale wavelet decomposition of conditional image generation) .

With multi-scale wavelet transformation, we can reformulate the conditional probability distribution of ${\bm{x}}_{0}$ on ${\bm{y}}$ as

$$ $p({\bm{x}}_{0}|{\bm{y}})=\Pi_{k=1}^{S}p({\bm{x}}_{H}^{k}|{\bm{x}}_{L}^{k},{\bm {y}}_{H}^{k})p({\bm{x}}_{L}^{S}|{\bm{y}}_{L}^{S}),$ (10) $$

where $S$ denotes the maximum scale and

$$ $({\bm{x}}^{1}_{H},{\bm{x}}^{1}_{L})^{T}=A{\bm{x}}_{0},\ ({\bm{x}}^{k+1}_{H},{ \bm{x}}^{k+1}_{L})^{T}=A{\bm{x}}^{k}_{L},\ k=1,\ldots$ (11) $$

Figure: Figure 2: Schematic diagram of MSCGM. The conditional image $y_{L}^{0}$ is first decomposed by multi-scale wavelet transform (WT). In the coarsest wavelet layer, a BBDP transforms the low-frequency subband of conditional image to the low-frequency subband of the target image. A multi-scale GAN transforms subsequent high-frequency subbands of conditional image to the high-frequency subbands of the target image and recovers the full-resolution image using inverse wavelet transform (IWT). $y_{L}^{i}$ and $y_{H}^{i}$ represent the low- and high-frequency wavelet coefficients of the conditional image at the $i^{th}$ level of wavelet transform, respectively. Similarly, $x_{L}^{i}$ and $x_{H}^{i}$ denote the low-and high-frequency wavelet coefficients of the target image at the $i^{th}$ level of wavelet transform.
Refer to caption: x2.png

Different from existing approaches, our method leverages BBDP and GANs to handle low- and high-frequency subbands at various scales respectively, and the schematic diagram of our model is illustrated in Fig. 2. For the coarsest level low-frequency subband ${\bm{x}}^{S}_{L}$, due to the whitening effect of the low-frequency subband after wavelet transform, DMs can effectively and efficiently approximate $p({\bm{x}}_{L}^{S}|{\bm{y}}_{L}^{S})$ with fewer sampling steps, while generating diverse and photorealistic images. For another, though the conditional distribution of high-frequency subbands deviates from unimodal Gaussian distributions considerably, the multi-scale GAN is able to approximate their multi-modal distribution and sample the full-resolution images rapidly in a coarse-to-fine style. Since the BBDP at the coarsest level produces samples with good diversity and fidelity, the possibility of mode collapse commonly observed in pure GAN models can be minimalized.

We adopt the Wasserstein distances between fake and real images [] to optimize the generator $G$ and discriminator $D$. The training loss for multi-scale GAN is

$$ $\mathscr{L}_{G}=\sum_{k=1}^{S}\left[\lambda(G({\bm{x}}_{L}^{k},{\bm{z}}^{k})-{ \bm{x}}_{H}^{k})^{2}+\nu(1-\mathrm{SSIM}(G({\bm{x}}_{L}^{k},{\bm{z}}^{k}),{\bm {x}}_{H}^{k}))-\alpha D(G({\bm{x}}_{L}^{k},{\bm{z}}^{k}))\right]$ (12) $$

$$ $\mathscr{L}_{D}=\sum_{k=1}^{S}\big{(}D(G({\bm{x}}_{L}^{k},{\bm{z}}^{k}))-D({ \bm{x}}_{H}^{k})\big{)}$ (13) $$

Here ${\bm{z}}^{k}$ refers to random white noise at scale $k$, $\mathrm{SSIM}(\cdot,\cdot)$ is the structural similarity index measure [].

The complete sampling process of our model is detailed in Appendix Algorithm 2.

###### Definition 4 (Multi-scale wavelet decomposition of conditional image generation) .

## 4 Experiments

In this section, we first elucidate the design and training details of our method, as well as the preparation of training and testing datasets. Then, we evaluate our method on various image restoration tasks in computational and microscopy imaging, and compare it with baseline methods. At the end of this section, we further perform ablation study to analyze the importance of essential components in MSCGM.

### 4.1 Experimental setup and implementation details

For the BBDP at the coarsest wavelet scale, we adapt the UNet architecture with multi-head attention layers as practiced in []. The number of sampling steps is set as 1000 for training. The Brownian bridge diffusion model (BBDM) baseline [] is implemented at the full resolution scale without wavelet decomposition, and the same 1000 discretization step is used for training. The training of IR-SDE (image restoration stochastic differential equation) [] and Refusion [] follow the training setups as their original setups.

Inspired by Chen et al. [], the generator adopt a similar architecture to NAFNet. The generator contains 36 NAFBlocks distributed at 4 scales. A $2\times 2$ convolutional layer with stride 2 doubling the channels connects adjacent scales in the encoding (downsampling) path, and a $1\times 1$ convolutional layer with pixel shuffle layer connects adjacent scales in the decoding (upsampling)superre path. The number of channels in the first NAFBlock is 64. The discriminator consists of 5 convolutional blocks and 2 dense layers at the end, and each convolutional block halves the spatial dimension but doubles the number of channels. The number of channels for the first convolutional block in the discriminator is 64. The details about training and dataset are elucidated in Appendix H.

### 4.2 Evaluation metrics

Peak Signal-to-Noise Ratio (PSNR) is commonly used to measure the quality of reconstruction in generated images, with higher values indicating better image quality. Structural Similarity Index Measure (SSIM) [] assesses the high-level quality of images by focusing on changes in structural information, luminance, and contrast. Fréchet Inception Distance (FID) score [] is used in generative models like GANs to compare the distribution of generated images against real ones, where lower FID values imply images more similar to real ones, indicating higher quality.

**Table 1: Comparison of IR-SDE, ReFusion, BBDM and our method (MSCGM) on $4\times$ super resolution experiment. Sample steps are set as 1000 for all methods. Metrics are calculated on $256\times 256$ center-cropped patches of DIV2K validation set, Set 5 and Set 14. Entries in bold indicate the best performance achieved among the compared methods. Sampling time is measured at the same resource cost.**
| Methods | Sampling Time (s)$\downarrow$ | PSNR (dB)$\uparrow$ | SSIM$\uparrow$ |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- |
| DIV2K | Set5 | Set14 | DIV2K | Set5 | Set14 |  |  |
| IR-SDE | 19.51 | 23.54 | 26.73 | 22.71 | 0.56 | 0.72 | 0.54 |
| ReFusion | 17.25 | 21.39 | 22.82 | 22.13 | 0.43 | 0.52 | 0.49 |
| BBDM | 32.29 | 31.50 | 31.60 | 30.39 | 0.66 | 0.76 | 0.68 |
| MSCGM | 2.28 | 31.66 | 32.33 | 30.79 | 0.72 | 0.85 | 0.71 |

### 4.3 Results and comparison

In this section, we first evaluate our method on two microscopy image restoration tasks with different samples and then on three natural image restoration tasks, encompassing various popular image restoration applications in computational imaging and computational microscopy. First, we apply our method to microscopy images, where the degradation process is complex and unknown. Given the pronounced contrast and sparsity inherent in microscopy images, it is crucial to use generative models capable of handling multi-modal distributions to adapt effectively to complex microscopy datasets. We utilize our method to perform super-resolution on microscopy images of nano-beads and HeLa cells, transforming diffraction-limited confocal images to achieve resolution beyond the optical diffraction limit and match the image quality of STED microscopy. Next, we assess the adaptability of our method to various image restorations tasks of natural images, including a $4\times$ super-resolution task on DIV2K dataset, a shadow removal task on natural images (ISTD dataset) and on a low-light image enhancement task on natural images (LOL dataset). Through comparison against competitive methods on various testbeds, we demonstrate the superior effectiveness and versatility of our method for image restoration.

Figure: Figure 3: Comparison of CMSR [] and our method on microscopy image super-resolution of nanobeads. (a) Input images captured by a confocal microscope, (b, c) SR outputs of CMSR and our method, and (d) ground truths captured by an STED microscope of the same FOV. (e-l) Zoom-in regions marked by the corresponding white boxes in (a-c). Cross-section intensity values along the dashed line are plotted.
Refer to caption: x3.png

Microscopy Image Super-resolution: We evaluate our method on microscopy image super-resolution tasks and compare it with existing generative models in this field. Unlike natural image super-resolution, the LR images are not downsampled but sampled at the same spatial frequency as the HR images. However, the LR images are limited by the optical diffraction limit, which is equivalent to a convolution operation on the HR images with a low-pass point spread function (PSF). We apply our method to confocal (LR) images of fluorescence nanobeads to evaluate its capability to overcome the optical diffraction limit (see Methods for sample and dataset details). Figure 3 illustrates one typical field-of-view (FOV) of nanobead samples ($\sim$ 20-nm) captured using confocal and STED microscopy. The dimensions of nanobeads are considerably smaller than the optical diffraction limit ($\sim$ 250-nm) and nearby beads cannot be distinguished in confocal images. MSCGM and a competitive cross-modality super-resolution (CMSR) model [] are trained on the same training data and learn to transform LR confocal images to match HR STED images. As illustrated in Fig. 3(a-d), our method effectively reconstructs super-resolved nanobeads beyond the optical diffraction limit and outperforms CMSR in terms of SSIM and PSNR. Moreover, as shown in Fig. 3(e-l), our method reconstructs the nanobeads with smaller diameters, better contrast, and more accurate intensity values, matching well with the ground truth STED images.

As another demonstration, we apply our method to fluorescence imaging of HeLa cells and present Fig. 4, Appendix Fig. 11, 12 to further support its success in microscopy image restoration.

Figure: Figure 4: Microscopy image super-resolution of our method on HeLa cells. (a) LR confocal image, (b) SR output image and (c) HR STED image of the same FOV.
Refer to caption: extracted/5715016/Figures/Hela_figure.png

Natural Image Super-resolution: We train the three methods on the DIV2K training dataset and then test them on the DIV2K validation set, Set5, and Set14. Table 1 quantifies the super-resolution performance in terms of PSNR and SSIM on the three test sets. The reported MSCGM scores indicate better restoration quality over competitive methods. Table 1 also presents the sampling time of each method under the same resource cost (GPU memory). Our method (MSCGM) achieves superior image generation results and exhibits a 10x improvement in average single-image generation time compared to IR-SDE, an 8x improvement compared to ReFusion, and a 16x improvement compared to BBDM. Additional visualization results of the three methods on the DIV2K validation set are displayed in Appendix Fig. 7, showing better reconstruction fidelity and perceptual quality.

**Table 2: Quantitative evaluation metrics of our method (MSCGM) and competitive methods on ISTD dataset. Metrics calculated on full-resolution images.**
| Methods | PSNR (dB)$\uparrow$ | SSIM$\uparrow$ |
| --- | --- | --- |
| DC-ShadowNet [] | 26.38 | 0.922 |
| ST-CGAN [] | 27.44 | 0.929 |
| DSC [] | 30.64 | 0.843 |
| DHAN[] | 27.88 | 0.921 |
| BMNet [] | 30.28 | 0.927 |
| ShadowFormer [] | 30.47 | 0.935 |
| BBDM [] | 30.54 | 0.910 |
| MSCGM(Ours) | 31.08 | 0.915 |

Moreover, we implement fast sampling with fewer sampling steps on the same super-resolution tasks, from 4 to 1000 steps, and depict the results in Appendix Fig. 8. Additional comparative samples, showcasing the performance of our model against BBDM at fewer sampling steps, are provided in Appendix Figs. 13 and 14.

Natural Image Shadow Removal: For image shadow removal task, we train our paradigm on the ISTD training dataset [] of 1330 image triplets (shadow image, mask and clean image) and evaluate it on the ISTD test set of 540 triplets. Table 2 summarizes the performance of our method against competitive methods, including DC-ShadowNet [], ST-CGAN [], DSC [], DHAN [], BMNet [], ShadowFormer [] and BBDM []. Appendix Fig. 10 further show the visualization results of our method and competitive methods. In summary, our method effectively restores high-quality images from shadowed and low-light images while minimizing the artifacts and inconsistency in the output images.

Low-light Natural Image Enhancement: In this task, we train and evaluate our paradigm on LOL dataset [], which contains 485 training and 15 testing image pairs. Figure 9 qualitatively showcases the results of our method in comparison with RetinexFormer [] and BBDM. Our method demonstrates better image enhancement performance with better color accuracy and less noisy details. This effectiveness of our method on image enhancement is further confirmed by the quantitative evaluation results against RetinexFormer, HWMNet [] and BBDM, as shown in Appendix Table 5.

Ablation Study: We first perform an ablation study to examine the effectiveness of the multi-scale generation of our method. For this purpose, we train a NAFNet-based GAN without the multi-scale wavelet decomposition on natural image super-resolution and shadow removal tasks and compare it with our method. It is important to note that the BBDP at the coarse level is inherent to the multi-scale generative mechanism and the NAFNet model performs image restoration in a single inference. As summarized in Table 3, our method achieves better reconstruction over the NAFNet model, confirming the effectiveness of the multi-scale generation and BBDP at the coarsest level, which provides our model with high sample diversity comparable to BBDM.

We then perform an ablation study to examine the effectiveness of the multi-scale generation of our method.

**Table 3: Abalation study among our method (MSCGM), BBDM and NAFNet-based GAN. Metrics were evaluated on DIV2K validation set and ISTD test set respectively. Sampling steps were set as 1000 for MSCGM and BBDM baseline.**
| Method | $4\times$ Super-resolution | Shadow Removal |  |  |
| --- | --- | --- | --- | --- |
|  | FID | PSNR | FID | PSNR |
| MSCGM | 148.33 | 31.66 | 47.19 | 31.08 |
| BBDM | 143.73 | 31.50 | 45.22 | 30.54 |
| NAFNet | 162.93 | 31.73 | 54.65 | 30.77 |

## 5 Conclusion

To address the limitations of exsiting diffusion models in conditional image restoration, we demonstrate a novel generative model for image restoration based on Brownian bridge process and multi-scale wavelet transform. By factorizing the image restoration process in the multi-scale wavelet domains, we utilize Brownian bridge diffusion process and generative adversarial networks to recover different wavelet subbands according to their distribution properties, consequently accelerate the sampling speed significantly and achieve high sample quality and diversity competitive to diffusion model baselines.

Limitations: In this work, while we mainly present multi-scale wavelet transformation using the common Harr bases, the future work could investigate designing more efficient wavelet or learned bases. Besides, the adversarial training employed in this study can be unstable and sensitive to training setups. In summary, our method provides a practical solution to issues with generative learning and facilitates the applications of advanced generative models on computational and microscopy imaging.

## Appendix A Code and data availability

The codes of our reported method is available at [https://anonymous.4open.science/r/MSCGM-E114](https://anonymous.4open.science/r/MSCGM-E114). The DIV2K dataset [] is obtained from [https://data.vision.ee.ethz.ch/cvl/DIV2K/](https://data.vision.ee.ethz.ch/cvl/DIV2K/), and the ISTD dataset [] is obtained from [https://github.com/DeepInsight-PCALab/ST-CGAN](https://github.com/DeepInsight-PCALab/ST-CGAN). The microscopy image datasets (nanobeads and HeLa cells) are requested from Wang et al. [] and partial demo images are uploaded to [https://anonymous.4open.science/r/MSCGM-E114](https://anonymous.4open.science/r/MSCGM-E114).

## Appendix B Score regularity for discretization

###### Theorem 5 .

Suppose the Gaussian distribution $p=\mathcal{N}(0,\Sigma)$ and distribution ${\tilde{p}}_{0}$ from time reversed SDE, the Kullback-Leibler divergence between $p$ and $p_{\tilde{0}}$ relates to the covariance matrix $\Sigma$ as:
$KL(p\parallel{\tilde{p}}_{0})\leq\Psi_{T}+\Psi_{\Delta t}+\Psi_{T,\Delta t}$, with:

$$ $\displaystyle\Psi_{T}$ $\displaystyle=f\left(e^{-4T}\left|\mathrm{Tr}\left((\Sigma-\mathrm{Id})\Sigma \right)\right|\right),$ (14) $\displaystyle\Psi_{\Delta t}$ $\displaystyle=f\left(\Delta t\left|\mathrm{Tr}\left(\Sigma^{-1}-\Sigma(\Sigma- \mathrm{Id})^{-1}\log(\Sigma)/2+(\mathrm{Id}-\Sigma^{-1})/3\right)\right| \right),$ (15) $\displaystyle\Psi_{T,{\Delta t}}$ $\displaystyle=o({\Delta t}+e^{-4T}),\qquad{\Delta t}\rightarrow 0,T\rightarrow+\infty$ (16) $$

where $f(t)=t-\log(1+t)$ and $d$ is the dimension of $\Sigma$, $\mathrm{Tr}\left(\Sigma\right)=d$.

###### Proposition 1 .

For any $\epsilon>0$, there exists $T,{\Delta t}\geq 0$ such that:

$$ $\displaystyle(1/d)(\Psi_{T}+\Psi_{\Delta t})\leq\epsilon,$ (17) $\displaystyle T/{\Delta t}\leq C\epsilon^{-2}\kappa^{3},$ (18) $$

where $C\geq 0$ is a universal constant, and $\kappa$ is the condition number of $\Sigma$.

[] provides the proof outline for Theorem 5, based on the following Theorem 6,

###### Theorem 6 .

Let $N\in\mathbb{N}$, $\Delta t>0$, and $T=N\Delta t$. Then, we have that $\bar{x}_{t}^{N}\sim\mathcal{N}(\hat{\mu}_{N},\Sigma^{\widehat{N}})$ with

$$ $\displaystyle\Sigma^{\widehat{N}}$ $\displaystyle=\Sigma+\exp(-4T)\Sigma^{\widehat{T}}+\Delta t\Psi^{\widehat{T}}+ (\Delta t)^{2}R^{\widehat{T},\Delta t},$ (19) $\displaystyle\hat{\mu}_{N}$ $\displaystyle=\mu+\exp(-2T)\hat{\mu}_{T}+\Delta te^{\widehat{T}}+\frac{(\Delta t )^{2}}{2}r^{T,\Delta t},$ (20) $$

where $\Sigma^{\widehat{T}},\Psi^{\widehat{T}},R^{T,\Delta t}\in\mathbb{R}^{d\times d}$, $\hat{\mu}_{T},e^{\widehat{T}},r^{T,\Delta t}\in\mathbb{R}^{d}$, and $\|R^{T,\Delta t}\|+\|r^{T,\Delta t}\|\leq R$, not dependent on $T\geq 0$ and $\Delta t>0$. We have that

$$ $\displaystyle\Sigma^{\widehat{T}}$ $\displaystyle=-(\Sigma-\mathrm{Id})(\Sigma\Sigma^{-1})^{2},$ (21) $\displaystyle\Psi^{\widehat{T}}$ $\displaystyle=\mathrm{Id}-\frac{1}{2}\Sigma^{2}(\Sigma-\mathrm{Id})^{-1}\log( \Sigma)+\exp(-2T)\Psi^{\widetilde{T}}.$ (22) $$

In addition, we have

$$ $\displaystyle\hat{\mu}_{T}$ $\displaystyle=-\Sigma^{-1}T\Sigma\mu,$ (23) $\displaystyle e^{\widehat{T}}$ $\displaystyle=\left\{-2\Sigma^{-1}-\frac{1}{4}\Sigma(\Sigma-\mathrm{Id})^{-1} \log(\Sigma)\right\}\mu+\exp(-2T)\widetilde{\mu}_{T},$ (24) $$

with $\Psi^{\widetilde{T}},\widetilde{\mu}_{T}$ bounded and not dependent on $T$.

###### Theorem 7 .

Suppose that $\nabla\log p_{t}(x)$ is $\varphi^{2}$ in both $t$ and $x$ such that:

$$ $\displaystyle\sup_{x,t}\left\|\nabla^{2}\log p_{t}(x)\right\|\leq K,\qquad \quad\left\|\partial_{t}\nabla\log p_{t}(x)\right\|\leq Me^{-\alpha t}\|x\|$ (25) $$

for some $K$, $M$, $\alpha>0$. Then, $\|p-\tilde{p}_{0}\|_{\text{TV}}\leq\Psi_{T}+\Psi_{\Delta t}+\Psi_{T,\Delta t}$, where:

$$ $\displaystyle\Psi_{T}$ $\displaystyle=\sqrt{2}e^{-T}\operatorname{KL}\left(p\parallel\mathcal{N}(0, \mathrm{Id})\right)^{1/2}$ (26) $\displaystyle\Psi_{\Delta t}$ $\displaystyle=6\sqrt{\Delta t}\left[1+\mathbb{E}_{p}\left(\|x\|^{4}\right)^{1/ 4}\right]\left[1+K+M\left(1+1/2\alpha\right)^{1/2}\right]$ (27) $\displaystyle\Psi_{T,\Delta t}$ $\displaystyle=o\left(\sqrt{\Delta t}+e^{-T}\right)\qquad{\Delta t}\rightarrow 0 ,T\rightarrow+\infty$ (28) $$

Theorem 7 generalizes Theorem 5 to non-Gaussian processes. Please refer to [] for the complete proof.

###### Theorem 5 .

###### Proposition 1 .

###### Theorem 6 .

###### Theorem 7 .

## Appendix C Characteristics of high and low frequency coefficients in the wavelet domain

### C.1 Gaussian tendency of low-frequency coefficients in higher scales

In an image, pixel intensities are represented as random variables, with adjacent pixels exhibiting correlation due to their spatial proximity. This correlation often follows a power-law decay:

$$ $C(d)=\frac{1}{(1+\alpha d)^{\beta}},$ (29) $$

where $C(d)$ is the correlation between pixels separated by distance $d$, and $\alpha$ and $\beta$ characterize the rate of decay.

The wavelet transform (i.e., Haar wavelet transform), particularly its down-sampling step, increases the effective distance $d$ among pixels, thereby reducing their original spatial correlation. This reduction is crucial for applying the generalized Central Limit Theorem [, , , ], which requires that the individual variables (pixels, in this case) are not strongly correlated.

At scale $k$ in the wavelet decomposition, the low-frequency coefficients, $\bar{X}_{k}$, representing the average intensity over $n_{k}$ pixels, are calculated as:

$$ $\bar{X}_{k}=\frac{1}{n_{k}}(X_{1}+X_{2}+\ldots+X_{n_{k}}),$ (30) $$

where $n_{k}$ is the number of pixels in each group at scale $k$.

As the scale increases, the effect of averaging over larger groups of pixels, combined with the reduced correlation due to down-sampling, leads to a scenario where the generalized Central Limit Theorem can be applied. Consequently, the distribution of $\bar{X}_{k}$ tends towards a Gaussian distribution:

$$ $\bar{X}_{k}\xrightarrow{\text{d}}\mathcal{N}(\mu_{k},\frac{\sigma_{k}^{2}}{n_{ k}}),$ (31) $$

where $\mu_{k}$ and $\sigma_{k}^{2}$ are the mean and variance of the averaged intensities at scale $k$, respectively. This Gaussian tendency becomes more pronounced at higher scales due to the combination of reduced pixel correlation and the averaging process.

Figure: Figure 5: KL divergence between the standard normal distribution and normalized sample distribution with respect to the wavelet scale. Images were sampled from DIV2K dataset by $32\times 32$ and $64\times 64$ patches.
Refer to caption: extracted/5715016/Figures/KL_div_scale.png

In Fig. 5, we quantify the Gaussianity of low-frequency wavelet subbands at different scales using the Kullback-Leibler (KL) divergence, which measures the distance between the standard normal distribution and normalized sample distribution of the low-frequency wavelet coefficients. We sampled random patches of different resolutions ($32\times 32$ and $64\times 64$) from DIV2K dataset to calculate KL divergence. With the increasing of the scale, KL divergence decreases, validating the Gaussian tendency of low-frequency subbands after multi-scale wavelet transforms. Figure 6 further validates this tendency by plotting the kurtosis of sample distribution of microscopy images of nanobeads with respect to the scale.

### C.2 Sparsity and non-Gaussianity of high-frequency coefficients

High-frequency coefficients, when analyzed through wavelet transform, exhibit a distinct property of sparsity, characterized by a majority of wavelet coefficients being near or at zero, with only a sparse representation of significant non-zero coefficients. This sparsity highlights the efficiency of wavelet transforms in encoding signal details and abrupt changes. Furthermore, these high-frequency components often deviate from Gaussian distributions, tending towards leptokurtic distributions [] with higher peaks and heavier tails. This non-Gaussian nature suggests a concentration of energy in fewer coefficients and is crucial in applications like signal denoising and compression, where recognizing and preserving these vital characteristics is paramount.

In the following proposition, we theoretically show that the conditional distribution of ${\bm{x}}_{H}^{k}$ on ${\bm{x}}_{L}^{k}$ exhibits highly non-Gaussian properties and yields sparse samples. For a given image ${\bm{x}}$ and threshold $t$, the sparsity of its high-frequency coefficients at $k$-scale is defined as:

$$ $s({\bm{x}}_{H}^{k})=\frac{\lVert\textbf{1}\{{\bm{x}}_{H}^{k}\leq t\}\rVert}{L^ {2}},\ k=1,2,\ldots$ (32) $$

Here $\lVert\cdot\rVert$ is the norm counting the number of 1s in the vector.
In this way, we could estimate the expected sparsity of the true marginal distribution $p({\bm{x}}_{H}^{k})$. Considering that the LL coefficients with approximate Gaussian distribution given the whitening effect of wavelet decomposition, we have the following proposition.

###### Proposition 2 .

For a sufficiently large $k$, if the expected sparsity of ${\bm{x}}_{H}^{k}$ has a lower bound $\alpha$

$$ $\mathbb{E}(s({\bm{x}}_{H}^{k}))\geq\alpha,$ (33) $$

where $\alpha\in[0,1]$. Then the conditional expected sparsity of ${\bm{x}}_{H}^{k}$ on ${\bm{x}}_{L}^{k}$ is bounded by

$$ $\mathbb{E}(s({\bm{x}}_{H}^{k})|{\bm{x}}_{L}^{k})\geq\alpha-\varepsilon,$ (34) $$

where $\varepsilon>0$ is a small positive number determined by $k$.

###### Proof.

According to Eq. 31, for a sufficiently large $k$ we could assume that

$$ $\int|p({\bm{x}}_{L}^{k})-f_{k}({\bm{x}}_{L}^{k})|d{\bm{x}}_{L}^{k}\leq\varepsilon,$ (35) $$

where $f_{k}({\bm{x}}_{L}^{k})$ is the PDF of standard Gaussian distribution.
Notice that

$$ $\displaystyle\mathbb{E}(s({\bm{x}}_{H}^{k}))$ $\displaystyle=\iint s({\bm{x}}_{H}^{k})p(s({\bm{x}}_{H}^{k})|{\bm{x}}_{L}^{k}) p({\bm{x}}_{L}^{k})d{\bm{x}}_{L}^{k}ds$ (36) $\displaystyle=\int\mathbb{E}(s({\bm{x}}_{H}^{k})|{\bm{x}}_{L}^{k})p({\bm{x}}_{ L}^{k})d{\bm{x}}_{L}^{k}\geq\alpha$ (37) $$

Since $s$ is a bounded function in $[0,1]$, $\mathbb{E}(s({\bm{x}}_{H}^{k})|{\bm{x}}_{L}^{k})$ has an uniform lower bound with respect to all ${\bm{x}}_{L}^{k}$, denoted as $\alpha^{\prime}$. In other words, there exists $\alpha^{\prime}\in[0,1]$ such that

$$ $\mathbb{E}(s({\bm{x}}_{H}^{k})|{\bm{x}}_{L}^{k})\geq\alpha^{\prime},\ \forall{ \bm{x}}_{L}^{k}$ (38) $$

We can get

$$ $\displaystyle\begin{aligned} &\int\mathbb{E}(s({\bm{x}}_{H}^{k})|{\bm{x}}_{L}^ {k})p({\bm{x}}_{L}^{k})d{\bm{x}}_{L}^{k}\\ =&\int\mathbb{E}(s({\bm{x}}_{H}^{k})|{\bm{x}}_{L}^{k})f_{k}({\bm{x}}_{L}^{k})d {\bm{x}}_{L}^{k}\\ &+\int\mathbb{E}(s({\bm{x}}_{H}^{k})|{\bm{x}}_{L}^{k})(p({\bm{x}}_{L}^{k})-f_{ k}({\bm{x}}_{L}^{k}))d{\bm{x}}_{L}^{k}\\ \geq&\ \alpha^{\prime}\end{aligned}$ (39) $$

Similarly, it is easy to see that $1$ is a trivial uniform upper bound for $\mathbb{E}(s({\bm{x}}_{H}^{k})|{\bm{x}}_{L}^{k})$. Thus,

$$ $\mathbb{E}(s({\bm{x}}_{H}^{k})|{\bm{x}}_{L}^{k})\geq\alpha^{\prime}=\int \mathbb{E}(s({\bm{x}}_{H}^{k})|{\bm{x}}_{L}^{k})p({\bm{x}}_{L}^{k})d{\bm{x}}_{ L}^{k}-\varepsilon\geq\alpha-\varepsilon.$ (40) $$

∎

###### Proposition 2 .

###### Proof.

### C.3 Quantifying Non-Gaussianity of datasets

Third- and fourth-order sample cumulants, i.e., skewness and kurtosis, to quantify the non-Gaussianity of certain sample distributions []. The non-Gaussianity of high-frequency subbands can be evidenced by the kurtosis plot with respect to wavelet scales in Fig. 6. The kurtosis of high-frequency subbands of microscopy images increases with the wavelet scales, showing the high non-Gaussianity of the distribution of high-frequency coefficients.

Figure: Figure 6: Kurtosis of (left) low-frequency and (right) high-frequency subband coefficients with respect to the wavelet scale. Metrics were calculated on microscopy images of nanobeads.
Refer to caption: x4.png

In the following subsections, we first introduce the definitions of the two metrics and then evaluate and compare the non-Gaussianity of DIV2K and microscopy nanobead datasets.

#### C.3.1 Skewness ( γ 1 subscript 𝛾 1 \gamma_{1} italic_γ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT )

Skewness is a measure of the asymmetry of the probability distribution of a real-valued random variable. It quantifies how much the distribution deviates from a normal distribution in terms of asymmetry. The skewness value can be positive, zero, negative, or undefined. In a perfectly symmetrical distribution, skewness is zero. Positive skewness indicates a distribution with an extended tail on the right side, while negative skewness shows an extended tail on the left side.

The mathematical formula for skewness is given by:

$$ $\gamma_{1}=E\left[\left(\frac{X-\mu}{\sigma}\right)^{3}\right]$ (41) $$

where $X$ is the random variable, $\mu$ is the mean of $X$, $\sigma$ is the standard deviation of $X$, and $E$ denotes the expected value.

The greater the absolute value of the skewness, the higher the degree of non-Gaussianity in the distribution.

#### C.3.2 Kurtosis ( β 2 subscript 𝛽 2 \beta_{2} italic_β start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT )

Kurtosis is a measure of the "tailedness" of the probability distribution of a real-valued random variable. It provides insights into the shape of the distribution’s tails and peak. High kurtosis in a data set suggests a distribution with heavy tails and a sharper peak (leptokurtic), while low kurtosis indicates a distribution with lighter tails and a more flattened peak (platykurtic). Kurtosis is often compared to the normal distribution, which has a kurtosis of 3 (excess kurtosis of 0).

The formula for kurtosis is:

$$ $\beta_{2}=E\left[\left(\frac{X-\mu}{\sigma}\right)^{4}\right]-3$ (42) $$

where the variables represent the same as in the skewness formula.

These statistical measures, skewness ($\gamma_{1}$) and kurtosis ($\beta_{2}$), are crucial for quantifying and analyzing the non-Gaussianity in image data. They provide valuable insights into the distribution characteristics of image pixel intensities, particularly in highlighting deviations from the normal distribution.

The higher the kurtosis, the greater the degree of non-Gaussianity in the distribution, indicating a distribution with heavier tails than a normal distribution.

#### C.3.3 Non-Gaussianity of datasets

Here we examine the non-Gaussianity of the distribution of DIV2K training dataset and the microscopy nanobead dataset used in this work. Table 4 summarizes their skewnesses and kurtoses. We observe that microscopy images tend to have larger absolute values of skewness and kurtosis, confirming their highly non-Gaussian distribution and the high condition number for microscopy image restoration problem. As a result, standard DMs with the assumption that the distribution of target images is close to normal does not hold for microscopy images and may lead to performance degradation and excessive sampling time.

**Table 4: Quantitative measure of non-Gaussianity of DIV2K and microscopy image datasets.**
| Dataset | Skewness | Kurtosis |
| --- | --- | --- |
| DIV2K | 0.2193 | -1.0373 |
| Microscopy nanobeads | 2.5852 | 13.3429 |

## Appendix D Wavelet transform

Wavelet transforms are derived from a single prototype function known as the ‘mother wavelet’. This function undergoes various scaling and translation processes to generate a family of wavelets. The general form of a wavelet function $\psi(x)$, derived from the mother wavelet, is expressed as:

$$ $\psi_{a,b}(x)=\frac{1}{\sqrt{\lvert a\rvert}}\psi\left(\frac{x-b}{a}\right),$ (43) $$

where $a,b\in\mathbb{R}$, $a\neq 0,\text{ and }x\in D$, with $a$ and $b$ representing scaling and translation parameters, respectively, and $D$ being the domain of the wavelets.

### D.1 Multiresolution analysis (MRA)

MRA is crucial in the Wavelet Transform [, ]. It involves constructing a basic functional basis in a subspace $V_{0}$ within $L^{2}(\mathbb{R})$ and then expanding this basis through scaling and translation to cover the entire space $L^{2}(\mathbb{R})$ for multiscale analysis.

For each $k\in\mathbb{Z}$ and $k\in\mathbb{N}$, we define the scale function space as ${\bf V}_{k}=\{f|f$ is restricted to $(2^{-k}l,2^{-k}(l+1))$ for all $l=0,1,\ldots,2^{k}-1$ and vanishes elsewhere$\}$. Each space ${\bf V}_{k}$ encompasses a dimension of $2^{k}$, with each subspace ${\bf V}_{i}$ nested within ${\bf V}_{i+1}$:

$$ ${\bf V}_{0}\subseteq{\bf V}_{1}\subseteq{\bf V}_{2}\subseteq\ldots\subseteq{ \bf V}_{n}\subseteq\ldots.$ (44) $$

Using a base function $\varphi(x)$ in ${\bf V}_{0}$, we can span ${\bf V}_{k}$ with $2^{n}$ functions derived from $\varphi(x)$ through scaling and translation:

$$ $\varphi_{l}^{k}(x)=2^{k/2}\varphi(2^{k}x-l),\quad l=0,1,\ldots,2^{k}-1.$ (45) $$

These functions, $\varphi_{l}(x)$, are known as scaling functions and they project any function into the approximation space ${\bf V}_{0}$. The orthogonal complement of ${\bf V}_{k}$ in ${\bf V}_{k+1}$ is the wavelet subspace ${\bf W}_{k}$, satisfying:

$$ ${\bf V}_{k}\oplus{\bf W}_{k}={\bf V}_{k+1},\quad{\bf V}_{k}\perp{\bf W}_{k}.$ (46) $$

Thus, we can decompose ${\bf V}_{k}$ as:

$$ ${\bf V}_{k}={\bf V}_{0}\oplus{\bf W}_{0}\oplus{\bf W}_{1}\oplus\ldots\oplus{ \bf W}_{k-1}.$ (47) $$

To form the orthogonal basis for ${\bf W}_{k}$, we construct it within $L^{2}(\mathbb{R})$. This basis is derived from the wavelet function $\psi(x)$, orthogonal to the scaling function $\varphi(x)$. The wavelet function is defined as:

$$ $\psi_{l}^{k}(x)=2^{k/2}\psi(2^{k}x-l),\quad l=0,1,\ldots,2^{k}-1.$ (48) $$

A key property of these wavelet functions is their orthogonality to polynomials of lower order, exemplified by the vanishing moments criterion for first-order polynomials:

$$ $\int_{-\infty}^{\infty}x\psi_{j}(x)dx=0.$ (49) $$

Orthogonality and vanishing moments are central to wavelets, enabling efficient data representation and feature extraction by capturing unique data characteristics without redundancy. This efficiency is particularly useful in areas like signal processing and image compression.

### D.2 Wavelet decomposition and reconstruction

The Discrete Wavelet Transform (DWT) provides a multi-resolution analysis of signals, useful in various applications [, , ]. For a discrete signal $f[n]$, DWT decomposes it into approximation coefficients $cA$ and detail coefficients $cD$ using the scaling function $\varphi(x)$ and the wavelet function $\psi(x)$, respectively:

$$ $cA[k]=\sum_{n}h[n-2k]f[n],$ (50) $$

$$ $cD[k]=\sum_{n}g[n-2k]f[n],$ (51) $$

where $h[n]$ and $g[n]$ are the low-pass and high-pass filters, respectively. This decomposition process can be recursively applied for deeper multi-level analysis.

Reconstruction of the signal $f^{\prime}[n]$ from its wavelet coefficients uses inverse transformations, employing synthesis filters $h_{0}[n]$ and $g_{0}[n]$:

$$ $\displaystyle f^{\prime}[n]=\sum_{k}cA[k]\cdot h_{0}[n-2k]+cD[k]\cdot g_{0}[n- 2k],$ (52) $$

where the synthesis filters are typically the time-reversed counterparts of the decomposition filters. In multi-level decompositions, reconstruction is a stepwise process, beginning with the coarsest approximation and progressively incorporating higher-level details until the original signal is reconstructed.

### D.3 Haar wavelet transform

The Haar wavelet [], known for its simplicity and orthogonality, is a fundamental tool in digital signal processing. Its straightforward nature makes it an ideal choice for a variety of applications, which is why we have incorporated it into our project. The discrete wavelet transform (DWT) using Haar wavelets allows for the efficient decomposition of an image into a coarse approximation of its main features and detailed components representing high-frequency aspects. This process, enhanced by multi-resolution analysis (MRA), facilitates the examination of the image at various scales, thereby uncovering more intricate details.

The mathematical representation of the Haar wavelet and its scaling function is as follows:

$$ $\displaystyle\psi(t)$ $\displaystyle=\begin{cases}1&\text{if }0\leq t<0.5,\\ -1&\text{if }0.5\leq t<1,\\ 0&\text{otherwise},\end{cases}$ (53) $$

$$ $\displaystyle\phi(t)$ $\displaystyle=\begin{cases}1&\text{if }0\leq t<1,\\ 0&\text{otherwise},\end{cases}$ (54) $$

where $\psi(t)$ denotes the Haar wavelet function and $\phi(t)$ is the corresponding scaling function.

The application of Haar wavelet transform extends to two-dimensional spaces, particularly in image processing. This extension utilizes the same decomposition approach as for one-dimensional signals but applies it to both rows and columns of the image. The filter coefficients for Haar wavelets are calculated through these inner product evaluations:

- •
Low-pass filter coefficients (h):
$\displaystyle h_{0}$
$\displaystyle=\int_{0}^{1}\phi(t)\cdot\phi(2t)\,dt,$
(55)
$\displaystyle h_{1}$
$\displaystyle=\int_{0}^{1}\phi(t)\cdot\phi(2t-1)\,dt.$
- •
High-pass filter coefficients (g):
$\displaystyle g_{0}$
$\displaystyle=\int_{0}^{1}\psi(t)\cdot\psi(2t)\,dt,$
(56)
$\displaystyle g_{1}$
$\displaystyle=-\int_{0}^{1}\psi(t)\cdot\psi(2t-1)\,dt.$

To normalize these coefficients (ensuring their L2 norm is 1), we find $h_{0}=\frac{1}{\sqrt{2}}$ and similar values for the other coefficients. Thus, the Haar filter coefficients are:

$$ $\displaystyle h$ $\displaystyle=\left[\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}\right],$ $\displaystyle\quad g$ $\displaystyle=\left[\frac{1}{\sqrt{2}},-\frac{1}{\sqrt{2}}\right].$ (57) $$

For two-dimensional DWT in image processing, these filter coefficients are matrix-operated on the image. The horizontal operation uses the outer product of the filter vector with a column vector, and the vertical operation uses the outer product with a row vector. The resulting filter matrices are:

$$ $\displaystyle H$ $\displaystyle=h^{T}\otimes h=\begin{bmatrix}\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{ 2}}\\ \frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\end{bmatrix},$ (58) $\displaystyle G$ $\displaystyle=g^{T}\otimes g=\begin{bmatrix}\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt {2}}\\ -\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\end{bmatrix}.$ $$

The matrix $H$ corresponds to low-pass filtering (approximation), while $G$ captures the high-frequency details. In image transformation using Haar wavelets, these matrices help derive various coefficients representing different aspects of the image, such as approximation($LL$), horizontal($LH$), vertical($HL$), and diagonal detail($HH$) components.

## Appendix E Duality proof

### E.1 Generative modeling in spatial domain

For the Score-based Generative Model (SGM) [, , ], the forward/noising process is mathematically formulated as the Ornstein-Uhlenbeck (OU) process. The general time-rescaled OU process is expressed as:

$$ $d{\bm{X}}_{t}=-g(t)^{2}{\bm{X}}_{t}dt+\sqrt{2}g(t)d{\bm{B}}_{t}.$ (59) $$

Here, ${({\bm{X}}_{t})}_{t\in[0,T]}$ represents the noising process starting with ${\bm{X}}_{0}$ sampled from the data distribution. ${({\bm{B}}_{t})}_{t\in[0,T]}$ denotes a standard d-dimensional Brownian motion. The reverse process, ${\bm{X}}^{\leftarrow}_{t}$, is defined such that ${({\bm{X}}^{\leftarrow}_{t})}_{t\in[0,T]}={({\bm{X}}_{T-t})}_{t\in[0,T]}$. Assuming $g(t)=1$ in standard diffusion models, the reverse process is:

$$ $d{\bm{X}}^{\leftarrow}_{t}=\left({\bm{X}}^{\leftarrow}_{t}+2\nabla\log p_{T-t} ({\bm{X}}^{\leftarrow}_{t})\right)dt+\sqrt{2}d{\bm{B}}_{t}.$ (60) $$

Here, $p_{t}$ is the marginal density of ${\bm{X}}_{t}$, and $\nabla\log p_{t}$ is the score. To revert ${\bm{X}}_{0}$ from ${\bm{X}}_{T}$ via the time-reversed SDE, accurate estimation of the score $\nabla\log p_{t}$ at each time $t$ is essential, alongside minimal error introduction during SDE discretization.

The reverse process approximation, as specified in Eq. 60, involves time discretization and score approximation $\nabla\log p_{t}$ by $s_{t}$, forming a Markov chain approximation of the time-reversed SDE. The chain starts with $\tilde{\mathbf{x}}_{T}\sim\mathcal{N}(0,I_{d})$, evolving over uniform time intervals $\Delta t$, from $t_{N}=T$ to $t_{0}=0$. The discretized process is detailed as:

$$ $\tilde{\bm{x}}_{t-1}=\tilde{\bm{x}}_{t}+\Delta t\left(\tilde{\bm{x}}_{t}+2s_{t }(\tilde{\bm{x}}_{t})\right)+\sqrt{2\Delta t}\mathbf{z}_{t},$ (61) $$

where $\mathbf{z}_{t}$ are instances of Brownian motion ${\bm{B}}_{t}$.

### E.2 Generative modeling in wavelet domain

Consider $\mathbf{X}$ as the image vector in the spatial domain. The discrete wavelet transform (DWT) [, ] of $\mathbf{X}$ can be expressed as:

$$ $\widehat{\mathbf{X}}=\mathbf{A}\mathbf{X},\quad\mathbf{X}\in\mathbb{R}^{d}.$ $$

Here, $\mathbf{A}$ represents the discrete wavelet transform matrix, which is orthogonal, satisfying $\mathbf{A}\mathbf{A}^{\top}=\mathbf{I}$. Various forms of $\mathbf{A}$ are utilized in practice, such as Haar wavelets.

In the context of score-based generative modeling, we consider the forward (or noising) process, which can be formulated by the Ornstein–Uhlenbeck (OU) process. This is mathematically described as:

$$ $d\mathbf{X}_{t}=-\gamma(t)^{2}\mathbf{X}_{t}dt+\sqrt{2}\gamma(t)d\mathbf{B}_{t},$ (62) $$

where $\mathbf{B}_{t}$ is a standard d-dimensional Brownian motion. Upon applying DWT to $\mathbf{X}_{t}$, the transformed $\widehat{\mathbf{X}}_{t}$ also follows the OU process:

$$ $d\widehat{\mathbf{X}}_{t}=-\gamma(t)^{2}\widehat{\mathbf{X}}_{t}dt+\sqrt{2} \gamma(t)\mathbf{A}d\mathbf{B}_{t},\quad\widehat{\mathbf{X}}_{0}=\mathbf{A} \mathbf{X}_{0}.$ (63) $$

Defining $\widehat{\mathbf{B}}_{t}=\mathbf{A}\mathbf{B}_{t}$, which also behaves as a standard Brownian motion. If $\mathbf{X}_{0}$ is sampled from a distribution $p$, then $\widehat{\mathbf{X}}_{0}$ originates from the distribution

$$ $\hat{q}=\mathcal{T}_{\mathbf{A}}\#p,$ (64) $$

where $\mathcal{T}_{\mathbf{A}}$ denotes the linear transformation operation by $\mathbf{A}$ and $\#$ represents the pushforward operation. Consequently, we have

$$ $\hat{q}(\mathbf{x})=p(\mathbf{A}^{\top}\mathbf{x}).$ (65) $$

Let $p_{t}$ be the density distribution of $\mathbf{X}_{t}$ and $\hat{q}_{t}$ that of $\widehat{\mathbf{X}}_{t}$. Then,

$$ $\hat{q}_{t}=\mathcal{T}_{\mathbf{A}}\#p_{t},\quad\hat{q}_{t}(\mathbf{x})=p_{t} (\mathbf{A}^{\top}\mathbf{x}).$ (66) $$

Define the score functions for both processes as:

$$ $\mathbf{s}_{t}=\nabla\log p_{t},\quad\mathbf{r}_{t}=\nabla\log\hat{q}_{t}.$ (67) $$

These functions are related by:

$$ $\mathbf{r}_{t}(\mathbf{x})=\frac{\nabla\hat{q}_{t}(\mathbf{x})}{\hat{q}_{t}( \mathbf{x})}=\frac{\mathbf{A}\nabla p_{t}(\mathbf{A}^{\top}\mathbf{x})}{p_{t}( \mathbf{A}^{\top}\mathbf{x})}=\mathbf{A}\mathbf{s}_{t}(\mathbf{A}^{\top} \mathbf{x}).$ (68) $$

For the reverse processes denoted as $\mathbf{X}^{\leftarrow}_{t}$ and $\widehat{\mathbf{X}}^{\leftarrow}_{t}$, assuming $\gamma(t)=1$, they are given by:

$$ $\displaystyle d\mathbf{X}^{\leftarrow}_{t}$ $\displaystyle=\left(\mathbf{X}^{\leftarrow}_{t}+2\mathbf{s}_{T-t}(\mathbf{X}^{ \leftarrow}_{t})\right)dt+\sqrt{2}d\mathbf{B}_{t},$ (69) $\displaystyle d\widehat{\mathbf{X}}^{\leftarrow}_{t}$ $\displaystyle=\left(\widehat{\mathbf{X}}^{\leftarrow}_{t}+2\mathbf{r}_{T-t}( \widehat{\mathbf{X}}^{\leftarrow}_{t})\right)dt+\sqrt{2}d\widehat{\mathbf{B}}_ {t}.$ $$

Here, $\widehat{\mathbf{B}}_{t}=\mathbf{A}\mathbf{B}_{t}$. Exploring the second SDE:

$$ $\displaystyle d\widehat{\mathbf{X}}^{\leftarrow}_{t}$ $\displaystyle=\left(\widehat{\mathbf{X}}^{\leftarrow}_{t}+2\mathbf{r}_{T-t}( \widehat{\mathbf{X}}^{\leftarrow}_{t})\right)dt+\sqrt{2}d\widehat{\mathbf{B}}_ {t}$ $\displaystyle=\left(\widehat{\mathbf{X}}^{\leftarrow}_{t}+2\mathbf{A}\mathbf{s }_{T-t}(\mathbf{A}^{\top}\widehat{\mathbf{X}}^{\leftarrow}_{t})\right)dt+\sqrt {2}d\widehat{\mathbf{B}}_{t}$ $\displaystyle\mathbf{A}^{\top}d\widehat{\mathbf{X}}^{\leftarrow}_{t}$ $\displaystyle=\left(\mathbf{A}^{\top}\widehat{\mathbf{X}}^{\leftarrow}_{t}+2 \mathbf{s}_{T-t}(\mathbf{A}^{\top}\widehat{\mathbf{X}}^{\leftarrow}_{t})\right )dt+\sqrt{2}\mathbf{A}^{\top}d\widehat{\mathbf{B}}_{t}.$ (70) $$

Substituting $\mathbf{A}^{\top}\widehat{\mathbf{X}}^{\leftarrow}_{t}$ with $\mathbf{X}^{\leftarrow}_{t}$ brings us back to the first equation. The training processes for $\mathbf{s}\theta,\mathbf{r}{\hat{\theta}}$ with $\mathbf{X}_{t}^{(i)},\widehat{\mathbf{X}}_{t}^{(i)}$ also follow the standard denoising score matching loss function:

$$ $\displaystyle\mathbb{E}_{t}\left\{\lambda(t)\mathbb{E}_{\mathbf{X}_{0}}\mathbb {E}_{\mathbf{X}_{t}|\mathbf{X}_{0}}\left[\left\|\mathbf{s}_{\theta}(\mathbf{X} _{t},t)-\nabla_{\mathbf{X}_{t}}\log p_{0t}(\mathbf{X}_{t}|\mathbf{X}_{0}) \right\|^{2}\right]\right\}$ $\displaystyle\mathbb{E}_{t}\left\{\hat{\lambda}(t)\mathbb{E}_{\widehat{\mathbf {X}}_{0}}\mathbb{E}_{\widehat{\mathbf{X}}_{t}|\widehat{\mathbf{X}}_{0}}\left[ \left\|\mathbf{r}_{\hat{\theta}}(\widehat{\mathbf{X}}_{t},t)-\nabla_{\widehat{ \mathbf{X}}_{t}}\log\hat{q}_{0t}(\widehat{\mathbf{X}}_{t}|\widehat{\mathbf{X}} _{0})\right\|^{2}\right]\right\}.$ (71) $$

The forward and reverse probability distribution functions $p_{0t}$ and $\hat{q}_{0t}$ are defined as per the standard SGM model:

$$ $\displaystyle p_{0t}(\mathbf{X}_{t}|\mathbf{X}_{0})$ $\displaystyle=\mathcal{N}(\mathbf{X}_{t};\sqrt{\overline{\alpha}_{t}}\mathbf{X }_{0},(1-\overline{\alpha}_{t})\mathbf{I}),$ $\displaystyle\hat{q}_{0t}(\widehat{\mathbf{X}}_{t}|\widehat{\mathbf{X}}_{0})$ $\displaystyle=\mathcal{N}(\widehat{\mathbf{X}}_{t};\sqrt{\overline{\alpha}_{t} }\widehat{\mathbf{X}}_{0},(1-\overline{\alpha}_{t})\mathbf{I}).$ (72) $$

## Appendix F Posterior distribution of Brownian bridge diffusion model

Brownian bridge process can prove to be Markovian and hence defined alternatively by its one-step forward process. The one-step forward process $q({\bm{x}}_{t}|{\bm{x}}_{t-1},{\bm{y}})$ can be derived as:

$$ $\displaystyle q({\bm{x}}_{t}|{\bm{x}}_{t-1},{\bm{y}})$ (73) $\displaystyle=$ $\displaystyle N(\frac{1-m_{t}}{1-m_{t-1}}{\bm{x}}_{t-1}+(m_{t}-\frac{1-m_{t}}{ 1-m_{t-1}}m_{t-1}){\bm{y}},\delta_{t|t-1}{\bm{I}}).$ $$

Here $\delta_{t|t-1}=\delta_{t}-\delta_{t-1}\frac{(1-m_{t})^{2}}{(1-m_{t-1})^{2}}$ and $m_{t}=\frac{t}{T}$.

The reverse process can be derived by Bayesian formula and shown to be a Gaussian process with mean $\mu^{\prime}_{t}(x_{t},x_{0},y)$ and variance $\delta_{t}^{\prime}{\bm{I}}$:

$$ $\displaystyle\begin{aligned} p({\bm{x}}_{t-1}|{\bm{x}}_{t},{\bm{x}}_{0},{\bm{y }})&=\frac{p({\bm{x}}_{t}|{\bm{x}}_{t-1},{\bm{y}})p({\bm{x}}_{t-1}|{\bm{x}}_{0 },{\bm{y}})}{p({\bm{x}}_{t}|{\bm{x}}_{0},{\bm{y}})}\\ &=N({\bm{x}}_{t-1};\bm{\mu}^{\prime}_{t}(x_{t},x_{0},y),\delta_{t}^{\prime}{ \bm{I}}),\end{aligned}$ (74) $\displaystyle\begin{aligned} \bm{\mu}^{\prime}_{t}({\bm{x}}_{t},{\bm{x}}_{0},{ \bm{y}})&=\frac{\delta_{t-1}}{\delta_{t}}\frac{1-m_{t}}{1-m_{t-1}}{\bm{x}}_{t} +(1-m_{t-1})\frac{\delta_{t|t-1}}{\delta_{t}}{\bm{x}}_{0}\\ &+(m_{t-1}-m_{t}\frac{\delta_{t-1}}{\delta_{t}}\frac{1-m_{t}}{1-m_{t-1}}){\bm{ y}},\end{aligned}$ (75) $\displaystyle\delta_{t}^{\prime}=\frac{\delta_{t|t-1}\delta_{t-1}}{\delta_{t}}.$ (76) $$

Utilizing the estimate of ${\bm{x}}_{0}$ by $\bm{\epsilon}_{\theta}$ in Eq. 5, we could eliminate ${\bm{x}}_{0}$ in Eq. 75 and rewrite Eq. 75 as []:

$$ $\displaystyle\begin{aligned} \bm{\mu}^{\prime}_{t}({\bm{x}}_{t},{\bm{y}})&= \frac{\delta_{t-1}}{\delta_{t}}\frac{1-m_{t}}{1-m_{t-1}}{\bm{x}}_{t}\\ &+(m_{t-1}-m_{t}\frac{\delta_{t-1}}{\delta_{t}}\frac{1-m_{t}}{1-m_{t-1}}){\bm{ y}}\\ &+(1-m_{t-1})\frac{\delta_{t|t-1}}{\delta_{t}}({\bm{x}}_{t}-\bm{\epsilon}_{ \theta}({\bm{x}}_{t},{\bm{y}},t))\\ &=c_{xt}{\bm{x}}_{t}+c_{yt}{\bm{y}}-c_{\epsilon t}\bm{\epsilon}_{\theta}({\bm{ x}}_{t},{\bm{y}},t),\end{aligned}$ (77) $\displaystyle\begin{aligned} c_{xt}&=\frac{\delta_{t-1}}{\delta_{t}}\frac{1-m_ {t}}{1-m_{t-1}}+(1-m_{t-1})\frac{\delta_{t|t-1}}{\delta_{t}}\\ c_{yt}&=m_{t-1}-m_{t}\frac{\delta_{t-1}}{\delta_{t}}\frac{1-m_{t}}{1-m_{t-1}} \\ c_{\epsilon t}&=(1-m_{t-1})\frac{\delta_{t|t-1}}{\delta_{t}}.\end{aligned}$ (78) $$

## Appendix G Training and sampling algorithms

Figure: Algorithm 1 Training of BBDP

Figure: Algorithm 2 Sampling

## Appendix H Datasets and implementation details

For super-resolution experiments on DIV2K, the LR images were bicubically downsampled to $64\times 64$ center-cropped patches and bilinearly upsampled to $256\times 256$ as the input images, and the HR images were center-cropped as $256\times 256$. For shadow removal on the ISTD dataset, we used a similar method, pairing original shadow images with their shadow-free counterparts, both consistently sized to maintain uniformity in processing.

For super-resolution experiments involving HeLa cell nuclei and nano-beads, we utilized a Leica TCS SP8 stimulated emission depletion (STED) confocal microscope, equipped with a Leica HC PL APO 100×/1.40-NA Oil STED White objective. The samples were prepared on high-performance coverslips (170$\mathrm{\SIUnitSymbolMicro m}$ $\pm$ 10$\mathrm{\SIUnitSymbolMicro m}$, Carl Zeiss Microscopy), facilitating precise imaging. In the staining process, HeLa cell nuclei were treated with Rabbit anti-Histone H3 and Atto-647N Goat anti-rabbit IgG antibodies. The nano-beads, used for the STED experiments, were processed with methanol and ProLong Diamond antifade reagents. Both samples were excited with a 633-nm wavelength laser. Emission detection was carried out using a HyD SMD photodetector (Leica Microsystems) through a 645–752-nm bandpass filter.

For super-resolution experiments on microscopy images of nanobeads, the low-resolution (LR) images were acquired with specific settings (16 times line average and 30 times frame average for nano-beads; 8 times line average and 6 times frame average for cell nuclei), ensuring the acquisition of high-quality images necessary for our computer vision analysis. The scanning step size, i.e., the effective pixel size is  30 nm. Here, all LR and HR images were center cropped as $256\times 256$ patches to match the setting in other experiments.

AdamW [] optimizers were used in all experiments with an initial learning rate of 1e-4. Models with exponential moving averaged parameters with a rate of 0.999 was saved and evaluated. The BBDM models at the coarsest wavelet scale were trained for 100000 steps with a batch size of 6.

In our GAN training, we used a batch size of 3, with images cropped at 256x256 resolution. We employed AdamW optimizers, setting the learning rate at 1e-4 for the generator and 1e-5 for the discriminator. The training loss function is a weighted sum of diverse loss terms: L1 loss with a weight of 20.0, adversarial loss at 0.1, and structural similarity index measure (SSIM) [] loss weighted at 0.5. The training was conducted for 200 epochs.

All models were implemented on NVIDIA V100 graphic cards using PyTorch. For speed test, a batch size of 16 was used for BBDM, IR-SDE and Refusion models, and a batch size of 64 was used for MSCGM so that they consume equivalent computational resources.

## Appendix I Additional results

### I.1 Visualization and fast sampling results on natural image superresolution

Figure: Figure 7: Comparison of our method and baselines on $4\times$ natural image super-resolution on DIV2K validation set. Sample steps are set as 1000 for all methods. Center crops of $256\times 256$ are shown.
Refer to caption: x5.png

Figure: Figure 8: Fast sampling comparison between our method (MSCGM) and BBDM on $4\times$ super-resolution. PSNR evaluated on center-cropped $256\times 256$ patches on DIV2K validation set.
Refer to caption: x6.png

### I.2 Results on image enhancement

Figure: Figure 9: Image enhancement results on LOL dataset. Center-cropped $256\times 256$ patches shown.
Refer to caption: x7.png

**Table 5: Comparison of HWMNet[], RetinexFormer[] and our method (MSCGM) on image enhancement experiment. Metrics calculated on $256\times 256$ center-cropped images of LOL dataset. Entries in bold indicate the best performance achieved among the compared method.**
| Methods | PSNR (dB)$\uparrow$ | SSIM$\uparrow$ | FID$\downarrow$ |
| --- | --- | --- | --- |
| BBDM | 28.66 | 0.860 | 98.49 |
| HWMNet | 29.40 | 0.902 | 82.30 |
| RetinexFormer | 29.55 | 0.897 | 109.75 |
| MSCGM | 29.27 | 0.863 | 115.96 |

### I.3 Visualization results on shadow removal task

Figure 10 visualizes the shadow removal performance of our approach and competitive methods. Overall, our method generates output images with better consistency and milder artifacts. Refer to Table 2 for quantitative assessments on these results.

Figure: Figure 10: Comparison of our method and baselines on the natural image shadow removal task on ISTD test set. Sample steps were set as 1000 for all methods.
Refer to caption: extracted/5715016/Figures/Shadow_figure.png

### I.4 Additional results on microscopy image super-resolution

Figure 4 and 11 presents additional results on microscopy image super-resolution of HeLa cells and Fig. 12 showcases additional results on microscopy image super-resolution of nanobeads. In correspondence to the results shown in the main text, SR images generated by our method match with the ground truths very well.

Figure: Figure 11: Additional microscopy image super-resolution results of HeLa cells.
Refer to caption: x8.png

Figure: Figure 12: Additional microscopy image super-resolution results of nanobeads.
Refer to caption: x9.png

### I.5 Fast sampling result

Here we present the sampling results of our method(MSCGM) and BBDM of each sampling steps (4, 16, 64, 256) on the Set14 and DIV2k validation dataset. Figure 13 show the sampling results of our method and BBDM with various sampling steps on Set14 dataset, and Fig. 14 illustrates the improvement of sampling quality with respect to the number of sampling steps from 4 to 1000 on DIV2K validation dataset. Remarkbly, our method recovers high-quality image considerably faster in fewer sampling steps, confirming its superiority in sampling speed compared to competitive diffusion models.

Figure: Figure 13: Sampling results of our method (MSCGM) and BBDM with various sampling steps from 4 to 1000. Images sampled from $64\times 64$ LR images in Set14.
Refer to caption: extracted/5715016/Figures/step4-1000_2.png

Figure: Figure 14: Sampling results of our method (MSCGM) and BBDM with various sampling steps from 4 to 1000. Images sampled from $64\times 64$ LR images in DIV2K validation dataset.
Refer to caption: extracted/5715016/Figures/step4-1000_1.png

## Peer Reviews

confidence: 5
contribution: 2
initial_rating: 6
presentation: 3
questions: 1. Could the authors provide more detailed explanations regarding the choice and role of each loss term in Equation 18 and explain how they determined the relative weighting (λ, ν, α values) between the terms.
2. Could the authors provide a comparison of training time and the number of training parameters for MSCGM versus other models?
3. Could the authors to provide a detailed algorithm or pseudocode for MSCGM training, similar to what they provided for BBDP Algorithm.
soundness: 3
strengths: 1. MSCGM’s wavelet-based decomposition and conditional modeling shows substantial improvements in sampling speed and better reconstruction quality.
2. By adapting the generative approach to frequency characteristics, MSCGM enhances detail in restored images, especially in high-frequency components crucial for microscopy images.
3. The authors presented a new loss function.
summary: The paper introduces a multi-scale conditional generative model (MSCGM) aimed at enhancing microscopic image restoration by combining wavelet transforms and a Brownian Bridge diffusion process. The authors leverage multi-scale wavelet transforms to efficiently model low- and high-frequency image components, significantly improving the generation quality and speed of image restoration compared to traditional diffusion models.
weaknesses: 1. Equation 18 combines multiple objectives—L2 loss, Structural Similarity Index Measure (SSIM), and Wasserstein distance—but the rationale behind each component’s inclusion is not fully explained. Additionally, the roles and relative importance of the scaling parameters λ, ν, and α are unclear. 
2. The training procedure for MSCGM is not explicitly described. Unlike the clear training steps outlined for BBDP, MSCGM lacks a step-by-step description of its training pipeline. 
3. While Table 1 compares MSCGM with other models in terms of PSNR, SSIM, and sampling time, it does not include training time or the number of trainable parameters for each method. Without these metrics, it is challenging to gauge MSCGM’s overall computational cost relative to other approaches. Including such details would provide a more comprehensive view of the model’s efficiency.
4. In Section 4.2, the authors state that FID is considered as an evaluation metric. However, this metric is not included in Table 1. As FID is widely used in assessing generative models for image quality, its inclusion would offer further insights into MSCGM’s performance in distributional similarity to real images.
5. Equations from 4 to 15 are borrowed from BBDP paper. It is better to include them under the Preliminaries section.

---

confidence: 5
contribution: 2
initial_rating: 3
presentation: 1
questions: Please see my concerns in Weakness
soundness: 2
strengths: 1. The authors recognize the loss of detail in LDM, a known issue, and apply it to the microscopic image restoration context, an interesting direction.
2. They introduce the novel idea that the Brownian bridge stochastic process could effectively integrate conditional images.
summary: The authors propose a multi-scale conditional generative model (MSCGM) for image restoration, incorporating multi-scale wavelet transforms and a Brownian bridge stochastic process. The wavelet transform is included due to its reversibility, which maintains information integrity in the latent diffusion space, in contrast to traditional Latent Diffusion Models (LDM). The Brownian bridge stochastic process is leveraged to introduce conditional images in both forward and reverse processes. While the authors aim to address microscopic image restoration, the motivation and results in the paper do not consistently support this focus.
weaknesses: 1. **Lack of Consistency:** The paper lacks organization and clarity. Although the title emphasizes "Microscopic Image Restoration," the experiments primarily focus on "Natural Image Super-resolution" and "Low-light Natural Image Enhancement." Only a small subset of results explores microscopic images. If the model is intended for general image restoration, it would be more accurate to propose it as a ‘unified image restoration’ model. I suggest the authors either refocus their experiments more heavily on microscopic image restoration to align with the title, or broaden the title to reflect the wider scope of image restoration tasks covered in the paper.
  
2. **Introduction Needs Refinement:** The introduction lacks a clear problem definition and research motivation. The first two paragraphs provide a broad overview of diffusion processes that diverges from the paper’s focus. The discussion on latent diffusion downsampling is a well-known issue and could be alleviated by higher resolutions. The authors should clearly articulate why microscopic images especially require the multi-scale wavelet transform in the introduction. Please include a discussion of how their approach compares to or builds upon these existing wavelet-based diffusion models in the Introduction, highlighting any key differences or improvements.
  
3. **Lack of Acknowledgment of Prior Work:** The paper does not credit previous studies applying wavelet transforms in diffusion models, which could mislead readers into believing the concept originated here. Papers like "Wavelet Diffusion Models are Fast and Scalable Image Generators (CVPR 2023)" and "Training Generative Image Super-Resolution Models by Wavelet-Domain Losses Enables Better Control of Artifacts (CVPR 2024)" are directly related and should be cited with comparisons to clarify this study’s contributions.
  
4. **Figure 1 Illustration Issues:** The paper title focuses on "Microscopic Image Restoration," yet Figure 1 uses natural images. Including examples of microscopic images to show the degradations introduced by LDM and Refusion compared to MSCGM would enhance clarity.
  
5. **Methodology Development Clarity:** The description of the wavelet transform on page 4 is overly general, with key details moved to the appendix. Clear explanations of any novel model designs or algorithmic adaptations should be provided in the main text.
  
6. **Quality of Mathematical Presentation:** Symbols in the equations are used without proper declarations or explanations. Inconsistent symbols, like the variable for the normal distribution \( N \), further detract from clarity.
  
7. **Algorithm 1 Lack of Context:** Algorithm 1 on page 5 is underdeveloped. Symbols are not defined before use, and the algorithm lacks defined input-output requirements.
  
8. **Figure 2 Diagram Confusion:** Figure 2 is difficult to interpret. The illustration doesn’t clearly label network modules, workflow processes, or shared parameters (only a line is shown), which fails to clarify the model structure effectively.
  
9. **Lack of Dataset Information:** The results section includes evaluations of microscopic images, but there’s no description of the dataset. Is it public or private? What is the image count? Without these details, readers cannot analyze or reproduce the results. Please provide a detailed description of the microscopic image dataset used, including its source, size, and any preprocessing steps applied.
  
10. **Insufficient Ablation Studies:** Results provide only a simple comparison with LDM, without deeper exploration of MSCGM’s components or ablation studies to justify the performance benefits of each module.
  
11. **Unconvincing Model Performance:** The model’s performance requires further validation through comparison with advanced models. Numerous diffusion-based image restoration models from 2024 exist, yet none are used for comparison. This weakens the paper’s credibility. Key diffusion-based image restoration works worth considering include:  
   - RDDM ([link](https://cvpr.thecvf.com/virtual/2024/poster/31373))  
   - HIR-Diff ([link](https://cvpr.thecvf.com/virtual/2024/poster/29665))  
   - WF-Diff ([link](https://cvpr.thecvf.com/virtual/2024/poster/30059))  
   - DeqIR ([link](https://cvpr.thecvf.com/virtual/2024/poster/31759))  
   - GDP ([link](https://cvpr.thecvf.com/virtual/2023/poster/22095))

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: Especially considering that inference time is one of the main benefits, why was it not compared to models with fewer step counts or at least an in-depth analysis of how step counts influence the SOTA model performance? E.g. [2], or other methods that can be applied to the problem domain?

[2] Phung, Hao, Quan Dao, and Anh Tran. "Wavelet diffusion models are fast and scalable image generators." Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. 2023.
soundness: 1
strengths: • Integrating the Brownian Bridge Diffusion Process (BBDP) with adversarial learning in a multi-scale wavelet diffusion framework is innovative, enhancing image quality and stability. 

• The model achieves notable speed improvements, delivering faster inference without sacrificing image quality. 

• Performance remains consistent across diverse experiments, demonstrating robustness on both microscopy and natural images.
summary: The authors present a novel multi-scale generative model that leverages the Brownian Bridge process within the wavelet domain. This approach enhances inference speed while maintaining high image quality and diversity. The model is further integrated with computational microscopy workflows, expanding its applicability. The authors evaluate its performance on both microscopy and natural image datasets, demonstrating that it achieves slightly better results compared to existing methods such as IR_SDE, Refusion, and BBDM, with the added advantage of faster inference.
weaknesses: • The paper lacks a clear motivation for applying this model to computational microscopy workflows. The rationale for this specific application is unclear and lacks context, the relevance to microscopy appears out of place. A discussion on how this functionality benefits microscopy would help justify this direction and clarify its practical utility.

• The primary advantage of this method is its reduced inference time; however, the paper lacks a direct comparison with other methods that similarly aim to improve efficiency. Including such a comparison would provide valuable context and help quantify the benefits more clearly.

• The general evaluation lacks depth and is missing ablation studies. 

• There appear to be configuration issues with the comparison methods. For instance, IR-SDE [1] is cited as requiring 100 steps, but the authors use 1000, which significantly prolongs inference time. With the correct configuration (100 steps), the inference time should drop from 32 seconds to approximately 3 seconds.

• The choice of metrics is limited and somewhat inadequate for a super-resolution task. Relying solely on PSNR and SSIM may overlook important aspects of image quality. Including pixel-based metrics would provide a more comprehensive evaluation and might show shortcomings of the proposed method.


[1] Luo, Ziwei, et al. "Image restoration with mean-reverting stochastic differential equations." arXiv preprint arXiv:2301.11699 (2023).

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: The paper demonstrated that the low-frequency coefficients in higher scales show Gaussian tendency and thus applied this to BBDP. The idea is novel and well hypothesized, but it would be helpful if other DM methods, such as IR-SDE and ReFusion methods that are implemented on 4x super-resolution experiment, are also tested on microscopy image dataset. Only CMSR (GAN: non-diffusion model), is compared at the moment, not showing the effectiveness of proposed near-Gaussianity assumption.
Similarly, applying BBDM to full resolution image does not seem to be fair comparison. Since many works demonstrated the effectiveness of multi-scale diffusion models, BBDM should be implemented in a same manner as the proposed method to prove the superiority of WT instead of other compression technique. Please conduct an ablation study that replaces WT with simple down-sampling.
Is there any specific reason why the proposed work adopted BBDM which was initially designed for image translation where input and target domains are different? Super-resolution tasks seem to have similar domains for input and target. Justify the choice of BBDM for super-resolution.
soundness: 3
strengths: The paper analyzed the characteristics of microscopic images and proposed adequate methodology to address the sparsity and non-Gaussianity. Since the wavelet transformation divides the image into two subbands (high- and low- frequency coefficients) losslessly, handling each subband in a different manner is original. 
The contribution of the work is clear and well demonstrated. 
In addition, the work could be further applied to different modality images where sparse or non-Gaussianity exist.
summary: The paper proposed a multi-scaled generative model that uses a diffusion model (DM) for low-frequency image and a GAN for high frequency images. The wavelet transform provides multi-scale images without lossy encoding process. The lossless compression is particularly important for microscopic imaging where high-frequency component are sparse and non-Gaussian. Additionally, the authors showed the near-Gaussian property of low-frequency component and thus employed Brownian Bridge Diffusion Process (BBDP). The idea of employing different networks (DM and GAN) to different resolutions according to the characteristics of microscopic dataset is novel. The proposed MSCGM (multi-scale conditional generative model) showed improved super-resolution result with fast inference time.
weaknesses: Although the idea of the paper is novel, the effectiveness of the work has not been thoroughly assessed. The use of WT and the superiority of the proposed method compared to conventional method should be further evaluated. The specific comments are described in Questions.

Meta-review: This work introduces a multi-scale generative model to enhance conditional image restoration by initiating the Brownian Bridge diffusion process specifically at the lowest-frequency subband and applying generative adversarial networks at subsequent multi-scale high-frequency subbands in the wavelet domain.

Decision: Reject