\section{Introduction}
Images captured under complicated rain scenarios inevitably suffer from the noticeable degradation of visual quality. The degradation causes detrimental impacts on many computer vision tasks such as image segmentation \cite{liu2020deep}, object detection \cite{wojna2019devil} and video surveillance \cite{sultani2018real}. Thus, it is indispensable to develop effective algorithms to recover high-quality rain-free images, which is referred to as \textit{image deraining}.

Image deraining is an ill-posed problem.
To make it well-posed, traditional methods usually exploit various image priors, such as sparse coding \cite{luo2015removing}, Gaussian mixture model \cite{li2016rain} and low-rank representation \cite{chang2017transformed}. However, these hand-crafted priors possess a limited representation ability,  leading to poor deraining results under complicated and varied rainy scenes.
By learning from massive synthetic clean/rainy image pairs, the performance of deep learning based techniques is substantially improved \cite{yang2017deep,fu2021rain}.
However, these supervised methods usually achieve sub-optimal performance on real-world rainy images, due to the fact that i) there exists the domain gap between synthetic and real-world rainy images, and ii) collection of real-world and large-scale clean/rainy image pairs is usually not possible.

To alleviate the aforementioned problems, semi-supervised deraining techniques, as promising solutions,
leverage paired synthetic data for good initialization and unpaired real-world data for generalization \cite{wei2019semi,ye2021closing,yasarla2020syn2real}. Unfortunately, their performance still depends on the domain gap between the synthetic and real-world rainy images. Once the gap is much larger, the deraining results would be most unsatisfactory.
Further, CycleGAN-based unsupervised methods are proposed for real-world image deraining, avoiding paired training data \cite{wei2021semi,han2020decomposed,jin2019unsupervised,zhu2019singe,chen2022unpaired}. But these methods are known as being difficult to train, due to their complex adversarial objectives. Thus, they easily suffer from a series of problems of premature convergence, model collapse, and robust overfitting, leading to image degradation including the loss of image details, remnant rain, halo artifacts, and/or color distortion.

Recently, various diffusion models \cite{sohl2015deep,ho2020denoising,song2019generative} have gained wide interest in computer vision. Compared to GANs, diffusion models provide a stable training process and are more powerful for modeling image pixel distribution.
However,
standard diffusion models usually focus on synthetic degradation, such as image inpainting, image colorization, and super-resolution, in which it is easy to simulate large-scale paired synthetic data rather than real-world ones.
No work to-date exploits diffusion models for real-world image deraining.

In this paper, we propose a novel unsupervised learning paradigm based on denoising diffusion probabilistic models (DDPM), called RainDiffusion, to tackle the unfavorable prevailing problem of real-world image deraining.
To enable unsupervised training for unpaired real-world data, RainDiffusion includes two main interactive branches:
Non-diffusive Translation Branch (NTB) and Diffusive Translation Branch (DTB).
To bypass the difficulty in unpaired training of standard diffusion models, NTB, only used in the model training phase, makes the best use of cycle-consistent architecture to generate initial clean/rainy image pairs for conditional diffusion model training.
Given initial image pairs, DTB designs two conditional diffusion modules for high-quality real-world image deraining and rain generation.
Instead of popular unsupervised adversarial networks, RainDiffusion exploits the powerful generative ability of DDPMs to mitigate this limitation of unstable training and significantly improve the quality of final derained images (see Fig. \ref{ab_reals}).
Experiments on several well-known benchmark datasets show that RainDiffusion outperforms existing un/semi-supervised methods and achieves comparable performance against fully-supervised ones.
The main contributions are as follows:
\begin{itemize}
\item \textit{First Diffusion-based Model for Real-world Image Deraining.} We propose a novel unsupervised image deraining method based on denoising diffusion probabilistic models, called RainDiffusion. RainDiffusion is a non-adversarial training paradigm, serving as a new standard bar for real-world image deraining.

\item \textit{Stable Non-GAN-based Unsupervised Training for Unpaired Real-world Data.}
RainDiffusion consists of two main interactive branches: Non-diffusive Translation Branch (NTB) and Diffusive Translation Branch (DTB). Through these two branches, RainDiffusion can remove the dependency on paired training data, and avoid training unstable adversarial objectives.

\item \textit{State-of-the-art Image Deraining Performance.} Extensive experiments show the superiority of RainDiffusion in real-world image deraining. Moreover, our approach performs favorably against existing semi-supervised, unsupervised, and even supervised methods on synthetic datasets.
\end{itemize}

\begin{figure*}[!t] \centering
	\includegraphics[width=1\linewidth]{pipline2.pdf}
	\caption{The pipeline of RainDiffusion. RainDiffusion consists of two main interactive branches: Non-diffusive Translation Branch (NTB) and Diffusive Translation Branch (DTB).
	In detail, NTB, only used in the model training phase, fully exploits cycle-consistent architecture to generate initial clean/rainy image pairs. Note that $G_{\phi}^{A}$ and $G_{\phi}^{B}$ are employed to produce the rainy and rain-free images, respectively. $x$, $y$, $x'$, $y'$, $x''$ and $y''$ refer to the
	given rain-free image, given rainy image, generated rainy image, generated rain-free image, reconstructed rain-free image, and reconstructed rainy image.
	Given initial image pairs, DTB leverages two conditional diffusion modules for translation between clean and rainy images. Especially, the deraining diffusion module learns a conditional reverse process $p_{\theta}^{A}(x_{0:T}|\Tilde{x})$ without modifying the diffusion process $q(x_{1:T}|x_{0})$ for $x$. The diffusion module for rain generation is similar to it, except it works oppositely.
    }
	\label{fig:network framework}
\end{figure*}
\section{Related Work}
\textbf{Single Image Deraining.} Early single image deraining methods employ several hand-crafted priors to restore the rainy images. For example, in \cite{luo2015removing,kang2011automatic}, it is assumed that rain streaks are sparse and of high frequency in the image, and thus the problem is translated into progressive image decomposition via dictionary learning. Li et al.~\cite{li2016rain} adopt the Gaussian Mixture Model (GMM) to accommodate multiple orientations and scales of rain streaks. In \cite{chen2013generalized,chang2017transformed}, the low-rank representations are incorporated into deraining. In \cite{gu2017joint} and \cite{wang2017hierarchical}, guided filter and analysis sparse representation are integrated into sparse coding for finer layer decomposition.

Recently, deep learning-based methods have been substantiated to be effective in image deraining \cite{fu2017clearing}. Fu et al.~\cite{fu2017clearing} first introduce an end-to-end residual CNN to simplify the learning processing. Network modules, such as dense block \cite{li2018non,wang2019erl}, recursive block \cite{fan2018residual,li2018recurrent} and dilated convolution \cite{deng2020detail}, and structures, such as RNN \cite{li2018recurrent,ren2019progressive}, GAN \cite{zhang2019image,chen2022unpaired} and multi-stream networks \cite{yang2017deep,zhang2018density}, are validated to be effective in image deraining.
Although we have witnessed promising deraining results on synthetic datasets, these learning-based efforts trained on synthetic images generalize poorly to real-world images, typically because of the obvious domain gap.

To bridge the domain gap between synthetic and real-world rainy images, Some semi-supervised frameworks have been proposed \cite{wei2019semi,ye2021closing,yasarla2020syn2real}. For example, Wei et al.~\cite{wei2019semi} first propose a semi-supervised rain removal framework where they adopt a likelihood term imposed on Gaussian Mixture Model (GMM) and minimize the Kullback-Leibler (KL) divergence between synthetic and real rain. Yasarla et al. \cite{yasarla2020syn2real} further propose a Gaussian Process-based semi-supervised learning framework to obtain better generalization performance. However, these supervised/semi-supervised methods still require paired data. Such paired data is very difficult or even impossible to collect under real-world rainy scenes.
Motivated by a popular image-to-image translation architecture, i.e., CycleGAN \cite{zhu2017unpaired}, recent works \cite{zhu2019singe,jin2019unsupervised,wei2021deraincyclegan,chen2022unpaired} attempt to exploit the improved CycleGAN architecture and constrained transfer learning to jointly the
rainy and rain-free image domains. However, the aforementioned methods all rely on their complex adversarial objective to develop their algorithms, which leads to the difficulty of stable training. Thus, they are also not capable of obtaining promising real-world image deraining results.

\textbf{Denoising Diffusion Probabilistic Models.} Recently, denoising diffusion probabilistic models (DDPM) \cite{sohl2015deep,ho2020denoising,song2019generative} have exhibited their powerful ability in many vision tasks such as image super-resolution \cite{saharia2022image}, text-to-image generation \cite{gu2022vector}, and image segmentation \cite{amit2021segdiff}. Different from other GANs, diffusion models have a stable training process and provide more diversity because they are likelihood-based.
However, regular diffusion models still require paired data for training, and are so far not considered for real-world image deraining. This motivates us to propose the first diffusion-based model for real-world image deraining.

\section{Proposed Method}
Fig. \ref{fig:network framework} presents the overall architecture of RainDiffusion. To the best of our knowledge, the proposed RainDiffusion is the first diffusion-based model for real-world image deraining. To enable unsupervised training for unpaired real-world data, it consists of two main interactive branches: Non-diffusive Translation Branch (NTB) and Diffusive Translation Branch (DTB).
Especially, NTB is used to guide the generation of initial clean/rainy image pairs.
Then, DTB inputs image pairs generated from NTB into two conditional diffusion modules for translation between clean and rainy images.
In the following, we first briefly introduce the denoising diffusion probabilistic models. After that, we elaborate on our NTB and DTB.

\subsection{Denoising Diffusion Probabilistic Models}
Denoising diffusion probabilistic models (DDPM) slowly corrupt the training data with Gaussian noise and learn to reverse the corruption as a generative model \cite{sohl2015deep,ho2020denoising,song2019generative}. In the forward process, Gaussian noise is added sequentially onto an input image $x_{0} \sim q(x_{0})$ over $T$ time steps according to the following Markovian process:
\begin{equation}
q(x_{t}|x_{t-1})= \mathcal{N}(x_{t};\sqrt{1-\beta_{t}}x_{t-1},\beta_{t}I), \forall t \in \{1,...,T\}
\end{equation}
\begin{equation}
q(x_{1:T}|x_{0})= \prod_{t=1}^{T} q(x_{t}|x_{t-1})
\end{equation}
where $\beta_{t}$ is noise variance, $\mathcal{N}(.)$ is a Gaussian distribution, $I$ is an identity covariance matrix with the same dimensions as the input image $x_{0}$. An important property of this recursive formulation is that it also allows the direct sampling of $x_{t}$, when $t$ is drawn from a uniform distribution,
i.e., $\forall t \sim \mathcal{U}(\{1,...,T\})$:
\begin{equation}
\begin{aligned}
q(x_{t}|x_{0})&= \mathcal{N}(x_{t};\sqrt{\overline{\alpha}_{t}}x_{0},(1-\overline{\alpha}_{t})I) \Longrightarrow \\
x_{t} &= \sqrt{\overline{\alpha}_{t}}x_{0}+\epsilon_{t} \sqrt{1-\overline{\alpha}_{t}}
\end{aligned}
\end{equation}
where $\epsilon_{t} \sim \mathcal{N}(0,I)$, $\alpha_{t}=1-\beta_{t}$ and $\overline{\alpha}_{t}=\prod_{i=1}^{t}\alpha_{i}$. If we have the original image $x_{0}$ and fix a variance schedule $\beta_{t}$, we can sample any noisy
version $x_{t}$ via a single step. Similarly, reverse diffusion also adopts a Markov chain from $x_{T}$ onto $x_{0}$, albeit each step aims to gradually denoise the samples. Under small $\beta_{t}$ and large $T$, the reverse transition probability between $x_{t}$ and
$x_{t-1}$ can be approximated as a Gaussian distribution:
\begin{equation}
\label{reverse}
p_{\theta}(x_{t-1}|x_{t})= \mathcal{N}(x_{t-1};\mu_{\theta}(x_{t},t),\sigma_{\theta}(x_{t},t))
\end{equation}
\begin{equation}
\label{reverse_0:T}
p_{\theta}(x_{0:T})= p(x_{T}) \prod_{t=1}^{T} p_{\theta}(x_{t-1}|x_{t})
\end{equation}
where the reverse process is parameterized by a neural network that estimates $\mu_{\theta}(x_{t},t)$ and $\sigma_{\theta}(x_{t},t)$. Especially, $\beta_{t}$ can be learned jointly with the model or just kept constant to ensure that $x_{t}$ approximately follows a standard normal distribution.

The model then can be trained by optimizing a variational bound on negative data log likelihood $\mathbb{E}_{q(x_{0})}[-\log p_{\theta}(x_{0})] \leq L_{\theta}$, which can be  performed by:
\begin{equation}
\begin{aligned}
L_{\theta}=\mathbb{E}_{q}[\begin{matrix}\underbrace{D_{KL}(q(x_{T}|x_{0})||p(x_{T}))}\\ L_{T}\end{matrix}\, \begin{matrix}\underbrace{-\log p_{\theta}(x_{0}|x_{1})}\\ L_{0}\end{matrix}\\
+\begin{matrix}\sum_{t>1}\underbrace{ D_{KL}(q(x_{t-1}|x_{t}),x_{0})||p_{\theta}(x_{t-1}|x_{t}))}\\ L_{t-1}\end{matrix}]
\end{aligned}
\end{equation}
where $D_{KL}$ denotes the Kullback-Leibler divergence between two probability distributions. A common parametrization focuses on $\mu_{\theta}(x_{t},t)$ while ignoring $\sigma_{\theta}(x_{t},t)$ \cite{ho2020denoising}, as follows:
\begin{equation}
\label{epsilon}
\mu_{\theta}(x_{t},t)= \frac{1}{\sqrt{\alpha_{t}}}(x_{t}-\frac{\beta_{t}}{\sqrt{1-\overline{\alpha}_{t}}}\epsilon_{\theta}(x_{t},t))
\end{equation}
In this case, the network can be used to estimate the added noise $\epsilon_{t}$ by minimizing the following loss:
\begin{equation}
\label{noise_loss}
L_{err}= \mathbb{E}_{x_{0},t,\epsilon_{t}\sim \mathcal{N}(0,I)}[||\epsilon_{t}- \epsilon_{\theta}(\sqrt{\overline{\alpha}_{t}}x_{0}+\epsilon_{t} \sqrt{1-\overline{\alpha}_{t}},t)||^{2}]
\end{equation}
During inference, reverse diffusion steps are performed starting from a random sample $x_{T} \sim \mathcal{N}(0,I)$. For each step $t \in \{T,...,1\}$, $\mu$ is derived using Eq. \ref{epsilon} based on the estimated $\epsilon_{\theta}$, and $x_{t-1}$ is sampled based on Eq. \ref{reverse}, which can be expressed in closed form:
\begin{equation}
\label{noise_loss2}
x_{t-1}=\frac{1}{\sqrt{\alpha_{t}}}(x_{t}-\frac{\beta_{t}}{\sqrt{1-\overline{\alpha}_{t}}}\epsilon_{\theta}(x_{t},t))+\sigma_{t}z
\end{equation}
where $z \sim \mathcal{N}(0,I)$), which resembles one step of sampling via Langevin dynamics \cite{welling2011bayesian}.

To enable high-quality image deraining, we should learn a conditional reverse process $p_{\theta}(x_{0:T}|\Tilde{x})$ without modifying the diffusion process $q(x_{1:T}|x_{0})$ for $x$, where $x_{0}$ and $\Tilde{x}$ represent a clean image and a rainy image. During the training phase, we sample $(x_{0},\Tilde{x}) \sim q(x_{0},\Tilde{x})$ from a paired data distribution and learn a conditional diffusion model. We input $\Tilde{x}$ to the reverse process as follows:
\begin{equation}
\label{reverse_rain}
p_{\theta}(x_{0:T}|\Tilde{x})= p(x_{T}) \prod_{t=1}^{T} p_{\theta}(x_{t-1}|x_{t},\Tilde{x})
\end{equation}
The noise estimator network in Eq. \ref{noise_loss}, \ref{noise_loss2} is also replaced by $\epsilon_{\theta}(x_{t},\Tilde{x},t)$.

\subsection{Non-diffusive Translation Branch (NTB)}
Although conditional diffusion models have shown state-of-the-art image-conditional data synthesis and restoration capabilities, the implementation of the diffusion models in Eq. \ref{reverse_rain} still requires paired clean/rainy images. We will tackle this problem by designing a new Non-diffusive Translation Branch (NTB) without requiring adversarial training. As shown in Fig. \ref{fig:network framework}, NTB introduces cycle-consistent architecture to generate initial clean/rainy image pairs for conditional diffusion model training. Since the training rainy image $y$ and clean image $x$ are unpaired, NTB employs two non-diffusive generators with parameters $\phi^{A,B}$ to obtain initial translation estimates:
\begin{equation}\centering
\begin{split}
x' = G_{\phi}^{A}(x)\qquad x'' = G_{\phi}^{B}(x')\\
y' = G_{\phi}^{B}(y)\qquad y'' = G_{\phi}^{A}(y')
\end{split}
\end{equation}
where $x$, $y$, $x'$, $y'$, $x''$ and $y''$ refer to the given rain-free image, given rainy image, generated rainy image, generated rain-free image, reconstructed rain-free image, and reconstructed rainy image. Notably, we adopt U-Net \cite{ronneberger2015u} as our non-diffusive generators. Then the cycle-consistency loss can be used to limit the space of generated samples and retain the contents of generated images:
\begin{equation}\label{eq:cyc}
\begin{aligned}
L_{cyc} = E_{x\sim P_{data(x)}}[\left \| x''-x \right \|_{1}]+E_{y\sim P_{data(y)}}[\left \| y''-y \right \|_{1}]
\end{aligned}
\end{equation}

\textbf{The difference between NTB and other circulatory structures.} In unsupervised learning, the circulatory structures with cycle-consistency loss are commonly used for model training \cite{chen2022unpaired,wei2021deraincyclegan,zhu2017unpaired}. The differences between NTB and others lie in the following two aspects: 1) NTB does not require any discriminators for adversarial training. This means that our NTB provides a more stable training process. 2) The generators of NTB are designed to generate paired clean/rainy images for diffusion model training. Thus, both non-diffusive generators $G_{\phi}^{A}$ and $G_{\phi}^{B}$ are only used in the training phase and do not involve the testing phase.

\begin{figure}[!t]
	\centering
	\includegraphics[width=1.0\linewidth]{derain_raing.pdf}
    \caption{Training RainDiffusion on unpaired data. From (a) to (d) are the input real-world rainy images, real-world clean images, generated derained images and generated rainy images, respectively.}
	\label{ab_real1}
\end{figure}

\begin{algorithm}[!t]
    \caption{RainDiffusion training}
    \label{alg:training}
    \begin{algorithmic}[1]
    \REQUIRE
    Unpaired clean image $x$ and rainy image $y$

        \REPEAT
        \STATE Randomly sample a binary patch mask $P_{i}$\\
        \STATE $x_{i}= Crop(P_{i}\circ x), y_{i}= Crop(P_{i}\circ y)$
        \STATE $x'_{i} = Crop(P_{i}\circ G_{\phi}^{A}(x)), x''_{i} = Crop(P_{i}\circ G_{\phi}^{B}(x'))$\\
        \STATE $y'_{i} = Crop(P_{i}\circ G_{\phi}^{B}(y)), y''_{i} = Crop(P_{i}\circ G_{\phi}^{A}(y'))$\\
        \STATE $t^{A}, t^{B} \sim$ Uniform $(\{ 1,...,T \})$\\
        \STATE $\epsilon^{A},\epsilon^{B} \sim \mathcal{N}(0,I)$\\
        \STATE Take gradient descent step on\\
        \STATE $\nabla_{\theta,\phi}[||\epsilon^{A}- \epsilon_{\theta}^{A}(\sqrt{\hat{\alpha}_{t^{A}}}x_{i}+\epsilon_{t^{A}} \sqrt{1-{\hat{\alpha}_{t^{A}}}},x'_{i},t^{A})||^{2}$\\
        \STATE $+||\epsilon^{B}- \epsilon_{\theta}^{B}(\sqrt{\hat{\alpha}_{t^{B}}}x'_{i}+\epsilon_{t^{B}} \sqrt{1-{\hat{\alpha}_{t^{B}}}},x''_{i},t^{B})||^{2}$\\
        \STATE $+||\epsilon^{B}- \epsilon_{\theta}^{B}(\sqrt{\hat{\alpha}_{t^{B}}}y_{i}+\epsilon_{t^{B}} \sqrt{1-{\hat{\alpha}_{t^{B}}}},y'_{i},t^{B})||^{2}$\\
        \STATE $+||\epsilon^{A}- \epsilon_{\theta}^{A}(\sqrt{\hat{\alpha}_{t^{A}}}y'_{i}+\epsilon_{t^{A}} \sqrt{1-{\hat{\alpha}_{t^{A}}}},y''_{i},t^{A})||^{2}$\\
        \STATE $+\lambda_{cyc}L_{cyc}]$
        \UNTIL{converged}\\
    \STATE \textbf{return} $\theta^{A},\theta^{B}$(optional)
    \end{algorithmic}
\end{algorithm}
\subsection{Diffusive Translation Branch (DTB)}
Following the advantage of denoising diffusion probabilistic models in the area of generative modeling, DTB also designs two conditional diffusion modules, to mitigate the limitation of unstable training and improve the generalization ability for real-world image deraining and rain generation (as shown in Fig. \ref{ab_real1}).

As shown in Fig. \ref{fig:network framework}, the deraining diffusion module learns a conditional reverse process $p_{\theta}^{A}(x_{0:T}|\Tilde{x})$ without modifying the diffusion process $q(x_{1:T}|x_{0})$ for $x$.
Since real-world rainy benchmarks usually consist of images with various sizes, we adopt the path-based DDPM \cite{ozdenizci2022restoring} as our backbone for size-agnostic image deraining. Notably, we also try other conditional diffusion models as our backbone, elaborated in Sec. \ref{Ablation_Study}.
Especially, we rewrite the conditional reverse process $p_{\theta}^{A}(x_{0:T}|\Tilde{x})$ as follows:
\begin{equation}
\label{crop_rain}
p_{\theta}(x_{0:T}^{i}|\Tilde{x}^{i})= p(x_{T}^{i}) \prod_{t=1}^{T} p_{\theta}(x_{t-1}^{i}|x_{t}^{i},|\Tilde{x}^{i})
\end{equation}
where $x_{0}^{i}= Crop(P_{i}\circ x_{0})$ and $\Tilde{x}^{i}= Crop(P_{i}\circ \Tilde{x})$ denote $p \times p$ patches from a training set image pair $(x_{0},\Tilde{x})$, and the $Crop(.)$ operation extracts the patch from the location indicated by $P_{i}$. During the training phase, we randomly sample (with uniform probability) the $p \times p$ patch location for $P_{i}$ within the complete range of image dimensions. Note that we use $p = 128$ pixel for $P_{i}$. The pipeline of the diffusion module for rain generation is similar to it, except it works oppositely. The training phase of our RainDiffusion is outlined in Algorithm \ref{alg:training}.

\begin{algorithm}[!t]
    \caption{RainDiffusion testing (image deraining)}
    \label{alg:test}
    \begin{algorithmic}[1]
    \REQUIRE
    Rainy image $\Tilde{x}$, conditional diffusion model $\epsilon_{\theta}^{A}(x_{t},\Tilde{x},t)$, number of implicit sampling steps $S$, dictionary of $D$ overlapping patch locations.
    \STATE $x_{t} \sim \mathcal{N}(0,I)$
    \FOR{each $i= S,...,1$}
        \STATE $t = (i-1)\cdot T/S+1$
        \STATE $t_{next} = (i-2)\cdot T/S+1$ \textbf{if} $i>1$ \textbf{else} 0
        \STATE $\widehat{\Omega}_{t} = 0, M=0$
        \FOR{each $d= 1,...,D$}
        \STATE $x_{t}^{d} = Crop(P_{d}\circ x_{t}), \Tilde{x}^{d}= Crop(P_{d}\circ \Tilde{x})$
        \STATE $\widehat{\Omega}_{t}=\widehat{\Omega}_{t}+P_{d}\cdot \epsilon_{\theta}(x_{t}^{d},\Tilde{x}^{d},t)$
        \STATE $M=M+P_{d}$
        \ENDFOR
        \STATE $\widehat{\Omega}_{t} = \widehat{\Omega}_{t}\oslash M$
        \STATE $x_{t} \Longleftarrow \sqrt{\overline{\alpha}_{t_{next}}}(\frac{x_{t}-\sqrt{1-\overline{\alpha}_{t}}\cdot\widehat{\Omega}_{t}}{\sqrt{\overline{\alpha}_{t}}})+\sqrt{1-\overline{\alpha}_{t_{next}}}\cdot \widehat{\Omega}_{t}$
    \ENDFOR
    \STATE\textbf{return} $x_{t}$
    \end{algorithmic}
\end{algorithm}
\begin{figure}[!t]
	\centering
	\includegraphics[width=1.0\linewidth]{testing.pdf}
    \caption{Illustration of testing phase of RainDiffusion detailed in Algorithm \ref{alg:test}. }
	\label{fig_test}
\end{figure}

\begin{table*}[!t]\footnotesize
\centering
\footnotesize
\caption{Comparisons on metrics (PSNR/SSIM) with the different types of deraining methods. \textcolor{red}{Red} and \textcolor{blue}{blue} colors are used to indicate the $1^{st}$ and $2^{nd}$ ranks, respectively. }
\label{synthetic_result}
\begin{tabular}{c|c|c|c|c|c|c|c}
\toprule
\multicolumn{2}{c|}{Datasets}                      & \multicolumn{2}{c|}{Rain200L} & \multicolumn{2}{c|}{Rain800} & \multicolumn{2}{c}{SPA-Data} \\
\midrule
\multicolumn{2}{c|}{Metrics}                       & PSNR          & SSIM         & PSNR         & SSIM         & PSNR          & SSIM         \\\midrule
\multirow{2}{*}{Prior-based}     & DSC \cite{luo2015removing}& 27.34         & 0.849        & 18.56        & 0.600        & 34.95         & 0.942        \\
                                 &GMM \cite{li2016rain}& 29.05         & 0.872        & 20.46        & 0.730        & 34.30         & 0.943        \\\midrule
\multirow{5}{*}{Supervised}      &DDN \cite{detail_layer}& 32.38         & 0.926        & 21.16        & 0.732        & 34.70         & 0.934        \\
                                 &DID-MDN \cite{zhang2018density}& \textcolor{blue}{35.40}         & 0.961        & 21.89        & 0.795        & 34.68         & 0.930        \\
                                 &SPA-Net \cite{wang2019spatial}& 31.59         & 0.919        & 24.37        & \textcolor{blue}{0.861}        & 35.13         & 0.944        \\
                                 &MSPFN \cite{jiang2020multi}& 32.40         & 0.933        & 25.52        & 0.830        & 35.21         & \textcolor{red}{0.951}        \\
                                 &AirNet \cite{li2022all}& 34.90 & \textcolor{blue}{0.969}         & 23.77& 0.833 & 33.45 & 0.924        \\
                                 \midrule
\multirow{3}{*}{Semi-supervised} &SIRR \cite{wei2019semi}& 34.47         & \textcolor{blue}{0.969 }       & 24.36        & 0.859        & 34.85         & 0.936        \\
                                 &Syn2Real \cite{yasarla2020syn2real}& 34.39         & 0.965        & 23.74        & 0.799        & 33.66         & 0.938        \\
                                 &JRGR \cite{ye2021closing}& 26.54         & 0.872        & 22.16        & 0.822        & 31.12         & 0.902        \\\midrule
\multirow{6}{*}{Unsupervised}    &CycleGAN \cite{zhu2017unpaired}& 24.61         & 0.834        & 23.95        & 0.819        & 22.40         & 0.860        \\
                                 &RR-GAN \cite{zhu2019singe}& /             & /            & 23.51        & 0.757        & /             & /            \\
                                 &DerainCycleGAN \cite{wei2021deraincyclegan}& 31.49         & 0.936        & 24.32        & 0.842        & 34.12         & 0.950        \\
                                 &DCD-GAN \cite{chen2022unpaired}& /             & /            & \textcolor{blue}{25.61}        & 0.813        & \textcolor{blue}{35.30}         & 0.943        \\
                                 &NLCL \cite{ye2022unsupervised}             &17.06             &0.653  &15.91            & 0.540             & 33.82         & 0.946        \\
                                 & \textbf{Ours}  & \textcolor{red}{37.02}              &\textcolor{red}{0.979}             &  \textcolor{red}{26.24}            &   \textcolor{red}{0.862}            & \textcolor{red}{36.96}              &\textcolor{blue}{0.948} \\
                                 \bottomrule
\end{tabular}
\end{table*}
The testing phase of RainDiffusion is illustrated in Fig. \ref{fig_test} and outlined in Algorithm \ref{alg:test}. We observe that a large $T$ inevitably leads to costly sampling, e.g., when $T = 1000$. To address this problem, we use an implicit sampling strategy \cite{song2020denoising} to accelerate our sampling process (lines 3-4 in Alg. \ref{alg:test}). Implicit sampling using a noise estimator network can be performed by
\begin{equation}
\label{implicit}
\begin{split}
x_{t-1}=&\sqrt{\overline{\alpha}_{t-1}}(\frac{x_{t}-\sqrt{1-\overline{\alpha}_{t}}\cdot \epsilon_{\theta}(x_{t},\Tilde{x},t)}{\sqrt{\overline{\alpha}_{t}}})\\
&+\sqrt{1-\overline{\alpha}_{t-1}}\cdot \epsilon_{\theta}(x_{t},\Tilde{x},t)
\end{split}
\end{equation}
during accelerated sampling we only needs a subsequence $\tau_{1},...,\tau_{S}$ of the complete ${1,...,T}$ timestep indices, which can be performed by
\begin{equation}
\label{implicit2}
\begin{split}
\tau_{i}=(i-1) \cdot T / S+1
\end{split}
\end{equation}
Given any denoising time step $t$, we first extract all overlapping $p \times p$ patches from a grid-like arranged parsing scheme to decompose the image $\Tilde{x}$ of arbitrary size, where $D$ represents the total number of extracted patches, defining a dictionary of overlapping patch locations (line 7 in Alg. \ref{alg:test}).
Then, we use $\epsilon_{\theta}(x_{t}^{d},\Tilde{x}^{d},t)$ to estimate the additive noise for all overlapping patch locations $d \in {1,...,D}$ and accumulate these overlapping noise estimates at their respective patch locations in a matrix $\widehat{\Omega}_{t}$ of same size as the whole image (line 8 in Alg. \ref{alg:test}).
Finally, we normalize $\widehat{\Omega}_{t}$ based on the number of received estimates for each pixel (line 11 in Alg. \ref{alg:test}), and perform an implicit sampling update using the smoothed whole-image noise estimate $\widehat{\Omega}_{t}$ (line 12 in Alg. \ref{alg:test}).

\section{Experiment}
We implement RainDiffusion using Pytorch 1.6 on a system with 11th Gen Intel(R) Core(TM) i7-11700F CPU and Nvidia GeForce RTX 3090 GPU. For optimizing RainDiffusion, we use the Adam optimizer with a min-batch size of 4 to train the paradigm, where the momentum parameters $ \beta_{1}$ and $\beta_{2}$ take the values of 0.5 and 0.999, respectively. Besides, the initial learning rate for both generators and discriminators is set to 0.00002. In the experiment, we set $\lambda_{cyc}$ to $1.0$. For training, a 128 $\times$128 image patch is randomly cropped from the original image (or its horizontal flipped version). During inference, only the noise estimator network responsible for the desired image deraining was used.
\subsection{Evaluation Settings and Datasets}
We qualitatively and quantitatively compare the performance of RainDiffusion with different types of deraining approaches, including supervised, semi-supervised and unsupervised methods. Supervised methods include DSC \cite{luo2015removing}, GMM \cite{li2016rain}, DDN \cite{detail_layer}, DID-MDN \cite{zhang2018density}, SPA-Net \cite{wang2019spatial}, MSPFN \cite{jiang2020multi} and AirNet \cite{li2022all}, while semi-supervised methods include SIRR \cite{wei2019semi}, Syn2Real \cite{yasarla2020syn2real}, and JRGR \cite{ye2021closing}. For unsupervised methods, CycleGAN \cite{zhu2017unpaired}, RR-GAN \cite{zhu2019singe}, DerainCycleGAN \cite{wei2021deraincyclegan}, DCD-GAN \cite{chen2022unpaired} and NLCL \cite{ye2022unsupervised} are adopted. As the code of the
methods \cite{zhu2019singe,wei2021deraincyclegan,chen2022unpaired} is not available, we refer to some comparison results presented in their paper. For other methods, we retrain them with the codes provided by authors.

As RainDiffusion is trained in an unsupervised learning manner, real-world clean and rainy images are utilized for training our model. For evaluation, both synthetic and real datasets are employed. The synthetic datasets are: 1) Rain200L \cite{yang2017deep} containing 1800 synthetic images for training and 200 images for testing; 2) Rain800 \cite{zhang2019image} containing 700 synthetic images for training and 192 images (100 synthetic images and 92 real-world images) for testing. The real datasets are: 1) SIRR-Data \cite{wei2019semi} containing 147 images without ground truths; 2) Google search with ``real-world rainy images" containing 300 images without ground truths; 3) SPA-Data \cite{wang2019spatial} containing 29500 high-quality rainy/clean image pairs (28500 images for training and 1000 images for testing). Finally, we adopt the average Peak Signal to Noise Ratio (PSNR) and Structural Similarity Index (SSIM) metrics for quantitative evaluation of the deraining images.

\begin{figure*}[!t]
	\centering
	\includegraphics[width=1\linewidth]{real_result.pdf}
    \caption{Comparison of deraining performance on real-world rainy images. It is clearly shown that most deraining methods fail to remove rain streaks in heavy rain. Comparatively, our method is more successful to remove rain accumulation and obtain the cleanest result with very clear details.}
	\label{morecompare}
\end{figure*}
\subsection{Comparison with State-of-the-arts}
\textbf{Comparison on Synthetic Datasets.}
Table \ref{synthetic_result} reports the quantitative results of 15 baseline models on synthetic datasets including Rain200L and Rain800. We have the following observations: 1) compared with two common types of unsupervised methods, i.e, CycleGAN-based \cite{zhu2017unpaired,zhu2019singe,wei2021deraincyclegan} and constractive learning \cite{ye2022unsupervised}, our method obtains higher values of PSNR and SSIM, which verifies the robustness of the proposed RainDiffusion. 2) There is an obvious gap between semi-supervised and supervised methods, even worse than unsupervised ones. 3) Our method can achieve competitive results to several existing supervised ones. More qualitative results on synthetic datasets are provided in the supplementary material.

\textbf{Comparison on Real-world Datasets.} We further compare our method with 15 baseline models on the SPA-Data dataset, which contains real-world rainy images with corresponding ground truth. As shown in Table~\ref{synthetic_result}, our method performs better than all semi-supervised and unsupervised methods and shows competitive advantages to the fully-supervised methods. Further, Fig. \ref{morecompare} illustrates the deraining performance of compared approaches on some real-world rainy images without ground truth, which are more complex and challenging. RainDiffusion successfully handles various rain streaks with different patterns, and achieves better visual results than all compared methods. Especially, the derained results of the unsupervised methods will cause color and structure distortion, while our method can better preserve the color and structure of the image. More qualitative results on real-world rainy images are provided in the supplementary material.

\begin{table}[!t]
\centering
\footnotesize
\caption{Ablation study for different loss functions. $L_{err}(x,x')$ indicates the noise estimator network in Eq. \ref{noise_loss} is replaced by $\epsilon_{\theta}(x_{t},x',t)$ (line 8 in Alg. \ref{alg:training}).}
\label{loss_ab}
\setlength{\tabcolsep}{1.5mm}{
\begin{tabular}{cccc}
\toprule
Model       & w/o $L_{err}(x,x')$& w/o $L_{err}(x',x'')$& w/o $L_{cyc}$\\
PSNR / SSIM &24.80/0.817     &25.36/0.842     &20.07/0.708      \\\midrule
Model       & w/o $L_{err}(y,y')$& w/o $L_{err}(y',y'')$& Ours \\
PSNR / SSIM &25.09/0.824     &25.78/0.851     & \textbf{26.24/0.872}   \\
\bottomrule
\end{tabular}
}
\end{table}

\begin{figure}[!t] \centering
	\includegraphics[width=0.8\linewidth]{loss_cyc.pdf}
	\caption{Effectiveness of the balanced weight $\lambda_{cyc}$.}
	\label{fig:weight_cyc}
\end{figure}

\begin{table}[!t]
\centering
\footnotesize
\caption{Quantitative comparisons with different diffusion-based models on the Rain800 dataset.}
\label{ab_ddpm}
\begin{tabular}{cccc}
\toprule
Method       & SR3&DDPM & WeatherDiffusion \\
PSNR / SSIM &26.12/0.866     &26.01/0.863     &\textbf{26.24/0.872}      \\
\bottomrule
\end{tabular}
\end{table}
\subsection{Ablation Study}
\label{Ablation_Study}

\textbf{The Effect of Loss Function.} We evaluate the effectiveness of our hybrid loss function on the Rain800 dataset. Especially, we remove one component to each configuration at one time. For fair comparision, the same training settings are kept for all models testing. As shown in Table \ref{loss_ab}, the full structure of RainDiffusion achieves the highest performance in terms of PSNR and SSIM, which indicates that all the components of RainDiffusion are beneficial for effective deraining. Moreover, we also evaluate the effectiveness of the balanced weight $\lambda_{cyc}$. As shown in Fig. \ref{fig:weight_cyc}, the balanced weight $\lambda_{cyc}$ is chosen experimentally, which optimizes our network effectively.

\textbf{The Effect of Different Diffusion Models.}
We also the effects of different diffusion models used in DTB, including improved DDPMs (WeatherDiffusion \cite{ozdenizci2022restoring} and SR3 \cite{saharia2022image}), and the standard DDPM \cite{ho2020denoising}. The above DDPMs are conditional diffusion models, which require paired data for training. As shown in Table \ref{ab_ddpm}, benefit from RainDiffusion, they both can achieve satisfactory performance on unsupervised image deraining.

\section{Conclusion}
In this paper, we do not plan to bridge the gap between the synthetic/real rainy images and depend on any adversarial objectives. We propose a new unsupervised learning paradigm based on denoising diffusion probabilistic models, dubbed RainDiffusion, to tackle the unfavorable prevailing problem of real-world image deraining. RainDiffusion consists of two main interactive branches: Non-diffusive Translation Branch (NTB) and Diffusive Translation Branch (DTB). Especially, NTB fully exploits cycle-consistent architecture to generate initial clean/rainy image pairs for conditional diffusion model training. Given initial image pairs, DTB leverages two conditional diffusion modules to progressively refine the desired output. Experiments on both synthetic datasets and real-world rainy images demonstrate that RainDiffusion performs favorably against current unsupervised and semi-supervised methods, and also possesses competitive advantages to the fully-supervised ones.

{\small
\bibliographystyle{ieee_fullname}


\section{}

\textit{Research Notes of the \href{https://aas.org}{American Astronomical Society}}
(\href{http://rnaas.aas.org}{RNAAS}) is a publication in the AAS portfolio
(alongside ApJ, AJ, ApJ Supplements, and ApJ Letters) through which authors can
promptly and briefly share materials of interest with the astronomical community
in a form that will be searchable via ADS and permanently archived.

The astronomical community has long faced a challenge in disseminating
information that may not meet the criteria for a traditional journal article.
There have generally been few options available for sharing works in progress,
comments and clarifications, null results, and timely reports of observations
(such as the spectrum of a supernova), as well as results that wouldn’t
traditionally merit a full paper (such as the discovery of a single exoplanet
or contributions to the monitoring of variable sources).

Launched in 2017, RNAAS was developed as a supported and long-term
communication channel for results such as these that would otherwise be
difficult to broadly disseminate to the professional community and persistently
archive for future reference.

Submissions to RNAAS should be brief communications - 1,000 words or fewer
\footnote{An easy way to count the number of words in a Research Note is to use
the \texttt{texcount} utility installed with most \latex\ installations. The
call  \texttt{texcount -incbib -v3 rnaas.tex}) gives 57 words in the front
matter and 493 words in the text/references/captions of this template. Another
option is by copying the words into MS/Word, and using ``Word Count'' under the
Tool tab.}, and no more than a single figure (e.g. Figure \ref{fig:1}) or table
(but not both) - and should be written in a style similar to that of a
traditional journal article, including references, where appropriate, but not
including an abstract.

Unlike the other journals in the AAS portfolio, RNAAS publications are not
peer reviewed; they are, however, reviewed by an editor for appropriateness
and format before publication. If accepted, RNAAS submissions are typically
published within 72 hours of manuscript receipt. Each RNAAS article is
issued a DOI and indexed by ADS \citep{2000A&AS..143...41K} to create a
long-term, citable record of work.

Articles can be submitted in \latex\ (preferably with the new "RNAAS"
style option in AASTeX v6.2), MS/Word, or via the direct submission in the
\href{http://www.authorea.com}{Authorea} or
\href{http://www.overleaf.com}{Overleaf} online collaborative editors.

Authors are expected to follow the AAS's ethics \citep{2006ApJ...652..847K},
including guidance on plagiarism \citep{2012AAS...21920404V}.

\begin{figure}[h!]
\begin{center}
\includegraphics[scale=0.85,angle=0]{aas.pdf}
\caption{Top page of the AAS Journals' website, \url{http://journals.aas.org},
on October 15, 2017.  Each RNAAS manuscript is only allowed one figure or
table (but not both). Including the
\href{http://journals.aas.org//authors/data.html\#DbF}{data behind the figure}
in a Note is encouraged, and the data will be provided as a link in the
published Note.\label{fig:1}}
\end{center}
\end{figure}

\acknowledgments

Acknowledge people, facilities, and software here but remember that this counts
against your 1000 word limit.

\section{Introduction}\label{sec:intro}

If a high-enough cosmological density of dark matter (DM) compact objects (COs) were to exist in the universe \citep{1991ComAp..15..139N}, detection of gravitational lensing events would be expected in the light curves of a single Gamma-Ray Burst (GRB) \citep{paczynski1987gravitational}. If a CO happens to lie near the line of sight of a bright source, then light from that source will be bent around it due to its gravitational effects. This phenomenon is known as gravitational lensing and can be used as a tool to probe CO DM. Examples of bright sources that have been used to study such an effect are quasars and galaxies. Sources that emit huge bursts of energy in a short amount of time such as GRBs fall under the category of transient sources. In the present work, GRBs are used as probes since they occur in abundance and can be found at cosmological distances.

When COs acting as gravitational lenses have masses in the range of $10^{6}$ $M_\odot$ to $10^{9}$ $M_\odot$, the gravitational lensing effect is known as millilensing \citep{nemiroff2001limits}. The term $'milli'$ is used when the angular separation of the lensed images is on the order of milli-arcseconds, even if the images are not resolved angularly. The CO masses that can be probed using the data of individual GRBs fall in the millilensing range. The lower end of this mass range is limited by the duration and the brightness of the GRBs whereas the upper limit is decided based on the available background data following each GRB.

One of the first papers to explore this mass range was by \citet{blaes1992using}. In this paper, we explore the possibility of detecting gravitational millilensing in the light curves of GRBs. The idea of detecting COs due to gravitational lensing of distant sources was first suggested by \citet{press1973method}. Detection of gravitational lenses from any temporally-resolved sources may provide an estimate of the fraction of the universe composed of COs within a specific mass range. Alternatively, a failure to detect them will impose constraints on the values of CO mass and cosmological density.

GRBs can be used as probes to search for gravitational lensing effects as they occur in abundance and can be found at cosmological distances.
The Burst and Transient Source Experiment (BATSE) \citep{fishman1989batse} onboard the Compton Gamma-Ray Observatory \citep{gehrels1993compton} has detected about 2700 GRBs and the currently ongoing Gamma-ray Burst Monitor (GBM) \citep{meegan2009fermi} on board the Fermi Gamma-ray Space Telescope \citep{michelson2010fermi} has detected about 3000 GRBs. Previous mililensing searches in BATSE GRBs were conducted by \citet{nemiroff2001limits}, and \citet{Marani1999GrbMillilensing}.
\citet{ougolnikov2001search} tried to find millilensing in 1512 BATSE GRBs and \citet{li2014search} also searched for lensing in BATSE GRBs. However, none of these searches resulted in the detection of any millilensing.

Recently, however, there have been several claims of gravitational millilensing in the light curves of single GRBs. The first claim published in the year 2021 was of GRB 950830 by \citet{paynter2021evidence}. Subsequently, two papers claiming millilensing in GRB 090717A were by \citet{kalantari2021imprints}, and \citet{2022ApJ...934..106K}. After that, two papers claiming millilensing in GRB 200716C were published by \citet{yang2021evidence}, and \citet{wang2021grb}.  All of these papers claim that the first episode of emission (hereafter called a pulse) near trigger time is gravitationally echoed by a second pulse immediately following the first pulse, thereby indicating that both pulses are gravitationally-lensed images of the same parent pulse. Even more recent claims of millilensing now exist, but they are not mentioned here because they are currently being analyzed.

There are certain criteria that need to be fulfilled to be sure that the two pulses are millilensed images of the same parent pulse and not an unrelated emission. These criteria are; 1) the emission should be significantly above the background, 2) the two pulses must have identical light curves within error margins, 3) Both the pulses should have identical spectra within error, and 4) both pulses must originate from the same location in the sky.

Following preliminary analyses by \citet{Mukherjee_2021}, \citet{mukherjee2021light}, and \citet{mukherjee2022light} on GRBs 950830, 090717A, and 200716C, we present here more detailed analyses of the gravitational millilensing hypotheses. In comparison to the research notes, in this paper, a modified version of the light curve similarity test as detailed in section \ref{subsec:ST} is proposed. Too many bins close to the background in some GRBs might give a low chi-square value--indicating similarity although the two pulse shapes close to the peak look very different and so these results cannot be trusted. This new method would allow comparing different regions of the two pulses in a GRB, keeping the brightness factor and the time delay between the arrival of the first and the second pulse the same. The method, therefore, helps to ensure the detection of a gravitational lens. Section \ref{subsec:Data} gives a detailed description of the data used. In Sections, \ref{subsec:ST} and \ref{subsec:HT} a light curve similarity test and a hardness similarity test are described respectively. Sections \ref{subsec:950830}, \ref{subsec:090717A}, and \ref{subsec:200716C} describe the results for GRB 950830, GRB 090717A, and GRB 200716C respectively. The conclusion and summary are given in section \ref{subsec:Conclusion}.

\section{Detectors and Data}\label{subsec:Data}
The GRBs analyzed here were detected by Fermi's GBM and Compton's BATSE. The BATSE detectors were housed in 8 modules located at the corners of the CGRO spacecraft. The primary detectors in each of the modules, the Large Area Detectors (LADs), are made up of NaI (Tl) crystals that operate at an energy range of 30 keV to 1000 keV. The GBM consists of 14 detector modules; out of these 12 are Sodium Iodide (NaI) detectors and two are Bismuth Germanate (BGO) detectors. The NaI detectors cover an energy range of 4 keV to 2000 keV, whereas the BGO detectors cover an energy range of 200 keV to 40 MeV.

The Time Tagged Events (TTE) BATSE data with a precision of 2 microseconds has been utilized here for analyzing BATSE-detected GRB 950830. The four BATSE energy channels utilized were channel 1 (20 - 60 keV), channel 2 (60 - 110 keV), channel 3 (110 - 320 keV), and channel 4 (320 - 2000 keV).

The 64-millisecond continuous time (CTIME) GBM data and the 2-microsecond TTE data containing individual detector counts were utilized for the GBM-detected GRB 090717A and GRB 200716C. The GBM CTIME data consists of counts divided into 8 different energy channels. The energy ranges for these channels are; channel 1 (4.5 - 11.8 keV), channel 2 (11.8 - 26.9 keV), channel 3 (26.9 - 50.4 keV), channel 4 (50.4 - 101.6 keV), channel 5 (101.6 - 293.8 keV), channel 6 (293.8 - 537.8 keV), channel 7 (537.8 - 983.3 keV), and channel 8 (983.3 - 2000 keV). The GBM TTE data given in energy/photon was divided into the same 8 energy channels as the CTIME data. However, only the channels that detected an obvious signal above the background were used in the analyses of GRB 090717A and GRB 200716C.

The energy bands already defined by BATSE and the Fermi GBM are used in the hardness similarity test described in section \ref{subsec:HT}. The selection of different energy bands might affect the hardness ratio values although not likely to modify the conclusions.

\section{Light Curve Similarity Test}\label{subsec:ST}
If two pulses in a GRB are gravitational lens images of a single parent pulse, their light curves should be statistically identical. A $\chi^2$ analysis is used here to test for light-curve similarity. To compute $\chi^2$, first, the values of two independent variables were determined; one is $r_{echo}$, the ratio in fluence between the two pulses, and the other one is $t_{offset}$, the time difference between the arrival of the first pulse and the arrival of the second pulse. Specifically, $r_{echo}$ was used to artificially decrease the fluence of the first pulse in order to match it to the fluence of the second pulse, and $t_{offset}$ was used to align the two pulses temporally. Following \citep{Nemiroff2000PulseStartConjecture} and \citep{Hakkila2009GrbPulseStart}, it was assumed that $t_{offset}$ is the same at all energies and at all times. Since gravitational lensing does not change the relative timings internal to images, all source images should have the same light curve to within an amplitude scale factor $r_{echo}$.

Given any $r_{echo}$, the reduced formula used for the $\chi^2$ is
\begin{equation}
    \chi^2 =  \sum_{t}\frac{(r_{echo} P_1(t) - P_2(t))^2}{r_{echo}^2 P_1(t) + P_2(t) + B_1(t) + B_2(t)},
\end{equation}
adapted from \citet{press1992numerical}, where
$P_1(t)$ is the number of counts in the first pulse over time,
$B_1(t)$ is the number of counts in the background underlying the first pulse over time,
$P_2(t)$ is the number of counts in the second pulse over time, and
$B_2(t)$ is the number of counts in the background underlying the second pulse over time.
The denominator represents the total variance for $P_1(t)$ and $P_2(t)$ and therefore incorporates the Poisson noise inherent in the data. Since in the numerator, $P_1(t)$ is reduced by $r_{echo}$, the variance term for $P_1(t)$ is $r_{echo}^2$ times $P_1(t)$.

The first step in our analysis procedure was to determine a background level for each GRB in each energy channel. This was done with a polynomial fit. Next, the counts from all the detectors and all energy channels that detected an obvious signal over the background were summed. A single background level was then found for the summed energies. Next, $t_{offset}$ and $r_{echo}$ were determined simultaneously by minimizing $\chi^2$ in a `Pulse Alignment' procedure.

To perform pulse alignment, first, the start and end times of each pulse were found separately. The start time of the first pulse was chosen to be the time when the summed counts in bins of a specific duration increased to over some sigma ($\sigma$) value above the background fit -- nearest the pulse peak. The pulse was then considered to continue until the summed counts dropped below the same $\sigma$ value. To be consistent with a possible gravitational lens interpretation, both pulses were taken to have the same duration. Then for a range of values of $r_{echo}$, the recorded counts following the first pulse were shifted in multiples of microseconds to find the $t_{offset}$ and $r_{echo}$ that minimized $\chi^2$.

Comparing the two pulses using a wide duration that encompasses most of each pulse may not be the best way to test the similarity between the two pulses. This is because when we compare long and faint pulse regions far from the peak (usually near the start and end of the pulses), it is likely that we will get a low value of $\chi^2$ (indicating higher chances of similarity) even though the central regions look very different. Therefore our primary light-curve similarity test involves a comparison using only bins near the pulse peak. To do this, first, the bin that contained the maximum number of counts ($peak_{bin}$) was determined for the first pulse. Then, several bins both before and after the $peak_{bin}$ were selected. This isolated a high-signal testable central region in the first pulse. Then the corresponding region of the second pulse was determined using $t_{offset}$ and the duration of the first pulse. Finally, the central regions for the two pulses were compared for light curve similarity using the $\chi^2$ formula mentioned above.

\section{Hardness Similarity Test} \label{subsec:HT}
The hardness test is based on the conjecture that the gravitational deflection of a photon does not change that photon's energy so that the gravitational lensing magnification of a source should be the same at every wavelength \citep{Paczynski1986GrbLensing}.

For this Hardness Similarity Test, we analyze most of the measured flux for each pulse. Specifically, we use the start and end times for the first pulse as determined by a designated $\sigma$ above background, as discussed in section~\ref{subsec:ST}. We then use the consistent start and end times for the second pulse as obtained from the Pulse Alignment procedure as discussed in section~\ref{subsec:ST}, thereby ensuring that $t_{offset}$ remains the same.

Given the background fits, pulse start times, and pulse durations, the counts above the background $P_{n,c}$ were determined, where the first subscript `n' refers to the pulse number and the second subscript `c' refers to the energy channel, in both pulses, and in all energy bands. The hardness ratio for each energy channel $r_c$ is found by taking the ratio of the summed counts in the second pulse $P_{2,c}$ to the summed counts in the first pulse $P_{1,c}$ and is given by $r_c = P_{2,c} / P_{1,c}$. When summing over all energy channels, $P_{n,all}$ is also determined, followed by $r_{all}$. The errors in all of these ratios were based on the Poisson noise inherent in the backgrounds as well as the pulses \citep{Sampling1977Cochran}.

\section{GRB 950830} \label{subsec:950830}

We analyzed the raw TTE data of BATSE which has a time resolution of 2 microseconds. This short GRB triggered three detectors: 5, 6, and 7, and is prominent in four different energy channels: 1, 2, 3, and 4. To maximize $\chi^2$ results, a time resolution of 0.004 seconds was chosen for all GRB 950830 analyses.

{\it Light Curve Similarity Test:}
Due to the short duration of GRB 950830, a flat background -- a zero-degree polynomial -- was found to accurately fit the background for each of the 4 energy channels. The counts and the background fit were summed across all utilized energies. The values of $t_{offset}$ and $r_{echo}$ that minimized $\chi^2$ were found to be 0.401 seconds and 0.711 respectively. The $peak_{bin}$ for pulse 1 was found to be -0.04 seconds relative to trigger time.

Visual inspection of the pulses showed that they differed the most nearest their peaks. We, therefore, tested several time windows around the pulse peaks to look for a maximum $\chi^2$ difference. Many time windows returned significant and similar $\chi^2$ values, but here we report on only one of them: the time window that involved 5-time bins before and 9-time bins after $peak_{bin}$. The duration of this time window was 0.06 seconds. The start time for pulse 1 was -0.060 seconds relative to the trigger time.

\begin{figure}
    \includegraphics[width=0.9 \columnwidth, angle=0]{950830_1.jpg}
    \caption{The light curve of GRB 950830 in 0.004-second time bins. ST and ET refer to pulse `Start Time' and `End  Time' respectively as shown by the vertical dotted lines.}
    \label{fig:950830_1}
\end{figure}

Figure~\ref{fig:950830_1} shows the regions that have been compared. A $\chi^2$ analysis found that the light curves of the two pulses differ at about the 3.0-$\sigma$ level (or 99.732 \%).

Comparing the pulses in different time resolutions but with a similar time window affects the $\chi^2$-derived $\sigma$ value. In general, smaller bin sizes -- which correspond to higher time resolution comparisons --  resulted in a lower $\sigma$ value because both pulses are closer to the background. In contrast, larger bin sizes -- which correspond to lower time resolution comparisons -- resulted in higher $\sigma$ values. The latter trend continued up until the bin size became close to the duration of the time comparison window when significant light-curve shape information was lost. The trend remains the same for all the other GRBs that follow. Due to this inherent property of binning, the results cannot be trusted in either of the extreme situations. Hence it is necessary to check the results with several different time bins, at least up to a resolution where the shape information is intact and yet the $\chi^2$ results are not deceiving.

{\it Hardness Similarity Test:} We first re-determined the start and end times of the first pulse from the summed counts across all energies. For this test, the start time was chosen to be the time when the summed counts in 0.004-s bins increased to over 1-$\sigma$ above the background fit -- nearest the pulse peak. The pulse was then considered to continue until the summed counts dropped to below 1-$\sigma$. The start time of Pulse 2 was found by adding $t_{offset}$ to the start time of Pulse 1. The end time of pulse 2 was chosen to be the time that maintained the same duration as pulse 1. The resulting starting times of the two pulses were determined to be -0.088 and 0.313 seconds relative to the BATSE trigger time, respectively. The duration of each pulse was found to be 0.184 seconds.

\begin{figure}
    \includegraphics[width=0.9 \columnwidth, angle=0]{950830_2.jpg}
    \caption{Count ratios between the two pulses as a function of energy channel for GRB 950830.}
    \label{fig:950830_2}
\end{figure}

As seen in Figure~\ref{fig:950830_2}, the $r_3$ value for channel 3 is somewhat lower than the other energy channels. A $\chi^2$ analysis found that the ratio $r_3$ between counts of the two pulses in BATSE energy channel 3 differs from $r_{all}$ at the 1.6-$\sigma$ level (89.289 \%).  Furthermore, channel 3 differs from channel 4 at the 1.8 $\sigma$ level (92.893 \%). From a $\chi^2$ test performed to check how different all the $r$ values were from $r_{all}$, it was found that the $r$ values differed from $r_{all}$ at about 1.9 $\sigma$ (or 93.92 \%).

\section{GRB 090717A} \label{subsec:090717A}
For GRB 090717A it was found that in Fermi GBM data, 5 of the 12 NaI detectors were pointed within 60 degrees of the burst and simultaneously detected a signal significantly above the background. These detectors were 0, 1, 2, 9, and 10. Fermi GBM energy channels 2, 3, 4, 5, and 6 were used because they recorded a prominent signal. To maximize $\chi^2$ results, a time resolution of 0.256 seconds was chosen for all analyses.

{\it Light Curve Similarity Test:} First, a third-degree polynomial was found to give an acceptable fit to the background of each of the 6 energy channels. Next, the counts and the background fit were summed across all energies. It was found that this background fit under pulse 2 did not change significantly, hence a constant background averaged over the pulse 2 region was assumed for simplicity. The $t_{offset}$ and $r_{echo}$ values that minimized $\chi^2$ were found to be 42.282 seconds and 0.542 respectively. The $peak_{bin}$ for pulse 1 was found to be 5.632 seconds. Almost all the time windows near the $peak_{bin}$ returned significant and similar $\chi^2$ values, but we will only report one of them here. The time window that involved 10 bins before and 17 bins after $peak_{bin}$ -- a total of 7.168 seconds -- gave optimal results. The resulting start time for pulse 1 was found to be 3.072 seconds relative to the trigger time.

\begin{figure}
    \includegraphics[width=0.9 \columnwidth, angle=0]{090717A_1.jpg}
    \caption{The light curve of GRB 090717A in 0.256-second time bins. ST and ET refer to pulse `Start Time' and `End  Time' respectively as shown by the vertical dotted lines.}
    \label{fig:090717A_1}
\end{figure}

Figure~\ref{fig:090717A_1} shows the regions that have been compared. The Light Curve Similarity test shows that the two pulses have dramatically different temporal shapes -- at above 8.3 $\sigma$ (or 99.999999999999989 \%).

{\it Hardness Similarity Test:} For this test, the start and end times of the first pulse was found from the summed counts across all energies. The start time was chosen to be the time when the summed counts in 0.256-s bins increased to over 1-$\sigma$ above the background fit -- nearest the pulse peak. The pulse was then considered to continue until the summed counts dropped to below 1-$\sigma$. The start time of Pulse 2 was found by adding $t_{offset}$ to the start time of Pulse 1. The end time of pulse 2 was chosen to be the time that maintained the same duration as pulse 1. The resulting start times of the two pulses were determined to be -1.024 and 41.258 seconds relative to the trigger time. The duration of each pulse was found to be 22.016 seconds.

\begin{figure}
    \includegraphics[width=0.9 \columnwidth, angle=0]{090717A_2.jpg}
    \caption{Count ratios between the two main pulses as a function of energy channel for GRB 090717A}
    \label{fig:090717A_2}
\end{figure}

Figure~\ref{fig:090717A_2} shows the $r$-values for the energy channels 2, 3, 4, 5, and 6, from left to right. As seen in this figure, the $r$ value for channels 2, 3, and 4, including the 1-$sigma$ error bars are different from the dotted line that depicts $r_{all}$. It was also found that the $r$ values differed from $r_{all}$ at about 1.7 $\sigma$ (or 90.953 \%).

Some specific energy channels differed from other specific energy channels at a relatively high formal significance. For example, it was found that $r_2$ differed from $r_3$ and $r_4$ separately at the 3.2 $\sigma$ level. As these discrepancies were found later from inspection of the data, they are not definitive in excluding the lensing hypothesis.

\section{GRB 200716C} \label{subsec:200716C}
For GRB 200716C it was found that in Fermi GBM data, 4 of the 12 NaI detectors were pointed within 60 degrees of the burst and simultaneously detected a signal significantly above the background. These detectors were 0, 1, 2, and 9. Energy channels 2, 3, 4, 5, and 6 were used for all the analyses as they recorded a prominent signal. The GRB was binned up to a time resolution of 0.008 seconds to maximize $\chi^2$.

{\it Light Curve Similarity Test:} First, a second-degree polynomial was used to fit the background of each of the 6 energy channels. Then the counts and the background fits were summed across all energies. Over the course of pulse 2, it was found that the background fit under pulse 2 did not vary significantly, hence a constant background averaged over the pulse 2 region was assumed for simplicity. The $t_{offset}$ and $r_{echo}$ values that minimized $\chi^2$ were found to be 1.909 seconds and 0.654 respectively. Next, the $peak_{bin}$ for pulse 1 was found to be 0.336 seconds relative to the trigger time. Almost all the time windows near the $peak_{bin}$ returned significant and similar $\chi^2$ values, but we will only report one of them here: the time window that involved 19 bins before and 5 bins after the $peak_{bin}$. The duration of this time window is 0.2 seconds. The pulse 1 start time was found to be 0.184 seconds relative to the trigger time.

\begin{figure}
    \includegraphics[width=0.9 \columnwidth, angle=0]{200716C_1.jpg}
    \caption{The light curve of GRB 200716C in 0.008-second time bins. ST and ET refer to pulse `Start Time' and `End  Time' respectively as shown by the vertical dotted lines.}
    \label{fig:200716C_1}
\end{figure}

Figure~\ref{fig:200716C_1} shows the regions that have been compared. The Light Curve similarity test shows that for a $r_{echo}$ of 0.654, the two pulses differ from each other at above 8.3 $\sigma$ confidence level (or 99.999999999999989 \%).

{\it Hardness Similarity Test:}  The start time for this test was chosen to be the time when the summed counts in 0.008-s bins increased to over 2-$\sigma$ above the background fit -- nearest the pulse peak. The pulse was then considered to continue until the summed counts dropped to below 2-$\sigma$. For this GRB, 2-$\sigma$ was used in place of the 1-$\sigma$ used for the other two GRBs because the first pulse did not sustain a return to within 1-$\sigma$ of background before the second pulse started. The start time of Pulse 2 was found by adding $t_{offset}$ to the start time of Pulse 1. The end time of pulse 2 was chosen to be the time that maintained the same duration as pulse 1. The resulting start times of the two pulses were determined to be 0.016 and 1.925 seconds relative to the trigger time. The duration of each pulse was found to be 0.56 seconds.

\begin{figure}
    \includegraphics[width=0.9 \columnwidth, angle=0]{200716C_2.jpg}
    \caption{Count ratios between the two main pulses as a function of energy channel for GRB 200716C}
    \label{fig:200716C_2}
\end{figure}

Figure~\ref{fig:200716C_2} shows the $r$-values for the energy channels 2, 3, 4, 5, and 6, from left to right. Although it was found that the $r_2$ value differs from $r_3$ by 1.7 $\sigma$, the $r$ values of other energy channels were all in rough agreement, not only with each other but with the combined $r_{all}$. We therefore conclude that the two pulses had hardnesses that were statistically consistent.

\section{Conclusion and Summary}\label{subsec:Conclusion}
In this paper, a search for gravitational millilensing in GRB 950830, GRB 090717A, and GRB 200716C was reported. With regard to the first GRB in this list, it was found that the light curves of the two pulses in GRB 950830 differ from each other at about the 3.0-$\sigma$ level (or 99.732 \%). It was also found that the hardness ratio or $r$ values differed from $r_{all}$ at about 1.9 $\sigma$ (or 93.92 \%).

The next burst of the three, GRB 090717A, was found to have light curves for the two bright pulses that were, at times, extremely different. Specifically, the light curves of the two pulses were found to be different at above 8.3-$\sigma$ level. In terms of hardness, the two main pulses were also found to be somewhat different: the $r$ values differed from $r_{all}$ at about 1.9 $\sigma$ (or 93.962 \%).

Last, our analysis of GRB 200716C showed that, again, the two bright pulses were extremely different in light curve shape. Specifically, the light curves of the two pulses were found to be different at above 8.3-$\sigma$ level. The hardnesses of the two main pulses in this GRB were found to be similar -- within an expected statistical range for each to be lensed images of the same parent pulse.

Hence we conclude that none of the three GRBs provide clear evidence in support of a gravitational lensing effect. In particular, the latter two bursts tested, GRB 090717A and GRB 200716C, even considering previously published analyses are very unlikely to involve multiple-image gravitational lensing.

We thank Michigan Technological University for their general support, and Joshua Wood for help with Fermi GBM Data Tools \citep{GbmDataTools} to find the detectors that pointed within 60 degrees of GRB 090717A and GRB 200716C separately. We sincerely thank the anonymous reviewers whose comments/suggestions helped improve the clarity of the paper.