# Radio-astronomical Image Reconstruction with Conditional Denoising Diffusion
Model

M. Drozdova  Department of Computer Science, University of Geneva, 7 route de
Drize, 1227 Carouge, Switzerland  V. Kinakh  Department of Computer Science,
University of Geneva, 7 route de Drize, 1227 Carouge, Switzerland  O. Bait
Geneva Observatory, University of Geneva, 51 Chemin Pegasi, 1290 Versoix,
Switzerland  O. Taran  Department of Computer Science, University of Geneva, 7
route de Drize, 1227 Carouge, Switzerland  E. Lastufka  Department of Computer
Science, University of Geneva, 7 route de Drize, 1227 Carouge, Switzerland  M.
Dessauges-Zavadsky  Geneva Observatory, University of Geneva, 51 Chemin
Pegasi, 1290 Versoix, Switzerland  T. Holotyak  Department of Computer
Science, University of Geneva, 7 route de Drize, 1227 Carouge, Switzerland  D.
Schaerer  Geneva Observatory, University of Geneva, 51 Chemin Pegasi, 1290
Versoix, Switzerland  S. Voloshynovskiy Corresponding authors:
<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Computer Science,
University of Geneva, 7 route de Drize, 1227 Carouge, Switzerland

###### Abstract

Reconstructing sky models from dirty radio images for accurate source
extraction, e.g. source localisation and flux estimation, is a complex yet
critical task. It has important applications in galaxy evolution studies at
high redshift, particularly in deep extragalactic fields using for example the
Atacama Large Millimetre Array (ALMA). With the development of large-scale
projects like the Square Kilometre Array (SKA), we anticipate the need for
more advanced source extraction methods. Existing techniques, such as CLEAN
and PyBDSF, currently struggle to effectively extract faint sources,
highlighting the necessity for the development of more precise and robust
methods. The success of source extraction task critically depends on the
quality and accuracy of image reconstruction. Since the imaging process
represents an information-lossy operator, the reconstruction is characterized
by uncertainty. The current study proposes the application of stochastic
neural networks for the direct reconstruction of sky models from dirty images.
This approach allows to localize radio sources and to determine their fluxes
with the corresponding uncertainties, providing a potential advancement in the
field of radio source characterization. We utilize a dataset of 10164 images
simulated with the CASA tool simalma based on the ALMA Cycle 5.3 antenna
configuration. We apply conditional Denoising Diffusion Probabilistic Models
(DDPMs) to directly reconstruct sky models from these dirty images, followed
by processing through Photutils to extract the coordinates and fluxes of the
sources. To test the robustness of the proposed model which was trained on a
fixed water vapour value, we examine its performance under varying levels of
water vapor. The proposed approach demonstrates state-of-the-art in source
localisation, achieving over 90% completeness at a signal-to-noise ratio (SNR)
as low as 2. Additionally, the described method offers an inherent measure of
prediction reliability thanks to the stochastic nature of the chosen model. In
terms of flux estimation, the proposed model excels past PyBDSF performance,
accurately extracting fluxes for 96% of the sources in the test set, a notable
improvement over CLEAN+ PyBDSF’s 57%. Conditional DDPMs is a powerful tool for
image-to-image translation, yielding accurate and robust characterisation of
radio sources, and outperforming existing methodologies. While this study
underscores its significant potential for applications in radio astronomy, we
also acknowledge certain limitations that accompany its usage, suggesting
directions for further refinement and research. The Python implementation of
the proposed framework is publicly available at
https://github.com/MariiaDrozdova/diffusion-for-sources-characterisation

##  1 Introduction

Figure 1: Simplified pipeline for real and simulated observations. Each pair
of antennas collects sky observations that yield UV-visibility data. These
data are then aggregated via a gridding process into an aggregated UV grid.
Finally, an inverse Fourier transform is applied along with the effect of a
primary beam to generate the dirty image.

The advent of the Square Kilometer Array (SKA) will revolutionize radio
astronomy at centimeter (cm) and millimeter (mm) wavelengths. Already, the
various SKA precursors, such as the Australian Square Kilometre Array
Pathfinder (ASKAP; Hotan et al., 2021), the Murchison Widefield Array (MWA;
Tingay et al., 2013), the LOw-Frequency ARray (LOFAR; van Haarlem et al.,
2013), and the MeerKAT radio telescope (Jonas and MeerKAT Team, 2016), have
started to produce interesting results. The Atacama Large Millimeter Array
(ALMA), owing to its excellent sensitivity, led to several interesting
studies. Some of the studies relevant to our current work are the deep
extragalactic continuum surveys, which focus on the gas and dust in high-
redshift galaxies; for example, the Reionization Era Bright Emission Line
Survey (REBELS; Bouwens et al., 2022), the ALMA SPECtroscopic Survey in the
Hubble Ultra-Deep Field (ASPECS; Walter et al., 2016), the ALMA Large Program
to Investigate [CII] at Early times (ALPINE Le Fèvre et al., 2020; Béthermin
et al., 2020; Faisst et al., 2020), and the GOODS-ALMA survey at 1.1 mm
(Franco et al., 2018) and the Automated Mining of the ALMA Archive in the
COSMOS Field (A3COSMOS) data
set111https://sites.google.com/view/a3cosmos/home?authuser=0 (Liu et al.,
2019).

These modern radio interferometers lie at the cutting edge of instrumentation
in terms of the receivers, correlators, low-noise amplifiers, high-speed data
links, and other technologies. They are also expected to produce data on
exabyte scales (Scaife, 2020), which will require novel and efficient data
reduction and analysis techniques. However, in order to convert the observed
visibilities to science-ready images, the processing of the resulting data
still relies on the traditional CLEAN algorithm (Högbom, 1974) for imaging and
deconvolution.

In radio astronomy, images are reconstructed using visibilities measured at
discrete locations depending on the antenna configuration and integration time
(aperture synthesis imaging). If the field of view of the antenna is
relatively small, these visibilities can be converted into “dirty images” by
taking an inverse Fourier transform. The dirty images are further improved to
remove the various artifacts coming from the side lobes using a deconvolution
process. This entire process of constructing the deconvolved images from the
observed visibilities is referred to as image reconstruction. The most popular
algorithm for such an imaging and deconvolution process is the CLEAN algorithm
(Högbom, 1974; Clark, 1980), which essentially relies on identifying intensity
peaks modeled as point sources in the dirty image. These peaks are iteratively
removed (or CLEANED) after convolving them with the dirty beam. In each
iteration, peaks of progressively lower intensity are CLEANED until a user-
defined threshold is achieved. Several variants of the CLEAN algorithm have
been introduced to mitigate its various limitations; for example, the
multifrequency synthesis algorithm (mfs; Conway et al., 1990) and multiterm
mfs (mtmfs; Rau and Cornwell, 2011) for wide-band imaging. These are widely
used by the community with radio interferometers, such as ALMA, Very Large
Array (VLA), MeerKAT, and so on. However, all CLEAN-based algorithms do not
address the fundamental issue of image reconstruction in radio astronomy. As
the visibilities are only measured at discrete locations, imaging in radio
astronomy is an ill-posed problem, such that the same visiblities can be
mapped to multiple sky models. Further, the CLEAN algorithm can benefit from
prior knowledge of source positions, which usually requires expert human
knowledge.

In addition to the CLEAN method, there exist other methods for reconstruction
of sky-model images from dirty images. Notably, Tikhonov regularization
(Tikhonov and Arsenin, 1977) and Wiener filtering (Wiener, 1949) represent two
classic deterministic methods. These iterative methods operate in an iterative
manner to minimize the fidelity loss with some priors on the expected
solution: $\mathcal{L}({\bf s})=||\bf z-H{\bf s}||^{2}_{2}+\lambda\Omega({\bf
s})$, where $\bf z$ is a dirty image, $\lambda$ is a regularization constant,
and $\Omega({\bf s})$ denotes the regularizer. Considering probabilistic
methods with a prior $p({\bf s})$, one can define $\Omega({\bf s})\propto-\ln
p({\bf s})$.

Another method used in radio astronomy is based on the maximum entropy method
(MEM; see e.g., Gull and Daniell 1978; Cornwell and Evans 1985; Narayan and
Nityananda 1986), whereby the distribution with the greatest entropy is chosen
among all the distributions that fit the data. Further techniques include bias
methods, the total variation prior (Rudin et al., 1992), and Gaussian mixture
methods (Bishop, 2007). These methods are often handcrafted and deterministic,
complementing the CLEAN method.

Various deep learning models have recently been proposed to solve the problem
of image reconstruction in radio astronomy. Gheller and Vazza (2022) use an
autoencoder network to “denoise” the dirty images, particularly for diffuse
extended emission. Connor et al. (2022) used a feed-forward neural network to
achieve super-resolution on already CLEANED images. Schmidt et al. (2022) used
a deep learning(DL)-based approach motivated by super-resolution models, where
they use a series of convolutional layers on the sampled amplitude and phases
of the visibilities to reconstruct the missing visibilities. Complementing
these deep learning models, Bluebild (Tolley et al., 2023) offers a novel
approach, leveraging functional principal component analysis (PCA) to
synthesize images by decomposing the sky into eigenimages with distinct energy
levels.

At the same time, unprecedented features of DL allow visibilities to be
directly translated into estimates of sky models or source locations. This
allows the user to skip the stage of dirty-image computation. Taran et al.
(2023) introduced a method for localizing radio sources from reduced UV
samples. The approach taken by these authors involves a two-stage pipeline:
initially, a fully connected network predicts a binary map where a value of
“1” signifies the presence of a source and “0” indicates the background.
Following this, an autoencoder refines the predictions. With the same dataset
as used in the present study, Taran et al. (2023) achieved a purity of 91% and
a completeness of 74%. A detailed comparison of our results and theirs can be
found in Section 4.

In the present work, we propose an alternative DL-based solution: we use a
denoising diffusion probabilistic model (DDPM; Sohl-Dickstein et al. 2015; Ho
et al. 2020) to reconstruct the sky model from dirty images and then localize
the sources. DDPMs are a family of generative models that use a diffusion
process to represent data distribution. In its basic unconditional
formulation, this type of model learns to map artificially generated noise
with the desired target distribution.

DDPMs have shown remarkable performance on a wide range of tasks:
unconditional image generation in Ho et al. (2020), conditional image
generation in Dhariwal and Nichol (2021), and text-to-image generation in
Ramesh et al. (2022). DDPMs are also applicable to other image-to-image
translation tasks, such as super-resolution, colorization, and inpainting
(Saharia et al. (2022)). Wang et al. (2023) demonstrate for the first time
that DDPM can be used for image reconstruction for radio astronomical images
trained on data generated for the Event Horizon Telescope (Akiyama et al.,
2019). In summary, the stochastic nature of DDPMs is beneficial in addressing
ill-posed problems where multiple solutions exist for the same observation.

The present paper is organized as follows. In Section 2, we provide a brief
overview of the simulated radio data used in the DDPM training. Section 3
details the proposed pipeline, encompassing preprocessing steps, the DDPM
model, post-processing techniques, and the range of metrics employed to
evaluate performance. Then, in Section 4, we present our findings in terms of
the performance of the DDPM approach in reconstruction, localization, and flux
estimation. Additionally, we analyze the performance in relation to S/N and
the varying noise amplitudes caused by water vapor. In Sections 5 and 6, we
consider the limitations of the approach and discuss potential directions for
future research.

##  2 Simulated data

We performed the entire training of the proposed DDPM model on simulated ALMA
data for the 12m array. These data are the same as those used in Taran et al.
(2023) for their DL model training. We refer to Taran et al. (2023) for
further details of the simulations. All the simulations are produced using the
Common Astronomy Software Applications (CASA) data processing software v6.2
(McMullin et al., 2007; Bean et al., 2022). Briefly, we chose a fixed ALMA
configuration, cycle 5.3. These simulations are designed to somewhat match the
real observations from the A3COSMOS dataset (Liu et al., 2019) with the aim
being to apply to them in the future. Therefore, all the simulated pointings
are randomly distributed around the COSMOS field (J2000 10h00m28.6s
+02d12m21.0s) within a radius of 1 deg. Each pointing corresponds to an
independent sky model. In the simulations, a sky model corresponds to an image
of 51.2$\times$51.2′′ across. The various sources in this sky model are
modeled as a 2D Gaussian overlaid on a zero background. The various source
properties, such as their sizes (major and minor axis of the Gaussian),
positions on the sky, and fluxes in each sky model are varied randomly. We
fixed the ALMA band to Band-6 (230 GHz) split into 240 channels with a channel
width of 7.8 MHz. We added the effect of noise from water vapor and receivers
using a precipitable water vapor (PWV) parameter value of 1.796. For each sky
model, we produced the corresponding visibilities for a total integration time
of 20 mins while the phasecenter of the sky model was under transit. We then
produced the corresponding dirty images for each pointing using the CASA task
tclean by setting the niter parameters equal to zero. This task does the
appropriate UV-gridding and inverse Fourier transform of the visibilities,
along with the effect of primary beam to produce the corresponding dirty
image. Figure 1 shows a block diagram of the simulation pipeline. In total, we
have 9164 simulated pointings (sky models) and the corresponding dirty images.
This comprises a total of 27632 sources, plus 1000 simulations without
sources.

Figure 2: Effect of $\gamma$ value on sky model preprocessing. The
preprocessing of the sky model via Equation 1 is significantly affected by the
$\gamma$ value. The proposed method modifies the intensities of the image
pixels to improve the network’s capacity for accurate reconstruction. By
employing a nonlinear root operation, nonzero values are pushed closer to one,
effectively making the histogram more uniform. However, this modification also
impacts the relative distances. Therefore, an optimal $\gamma$ value is
required to ensure a compromise between these two tendencies. The histogram
y-axis limits are set so that the height of the third bin comprises 80% of the
image.

In this study, these dirty images do the conditioning for the DDPM model,
which predicts the corresponding sky model. In order to study the effect of
varying noise, we also add simulations with the same configuration, but change
the PWV values from 0.472 to 5.186. We note that these additional simulations
are not used for training the model, but simply for testing the effect of
varying noise, as in the analysis of Taran et al. (2023).

##  3 Proposed method

###  3.1 Preprocessing

In this work, we apply machine learning based on DDPM to dirty images in order
to predict sky models. Both dirty images and sky models require preprocessing.
In machine learning, natural images are usually mapped to the $[-1,1]$ range
by a simple scaling. As we work with measurement data, which have their own
properties and statistics, we need to adjust the data-preprocessing step.

The sky models are non-negative and are typically very sparse and of very
small magnitude, of order $10^{-5}$ Jy. We would like to not only change the
range of values for sky models but make the distribution of the values
flatter. A sky model is what the model predicts, and so the loss functions
will be applied directly to it. In this study, we use mean square error (MSE)
due to the model choice. MSE loss works better with dense images (the ones
whose values have flatter histograms). We therefore need a transform that
makes the image histograms flatter and is also invertible, because after
getting predictions, we need to return to the original ranges.

We take inspiration from the preprocessing techniques used in particle physics
Finke et al. (2021) to face the sparsity and low signal magnitude. These use
blur and root functions. As the sky models already contain sources modeled as
Gaussians, we retain only the root function. We use the following scaling in
the experiments:

${\bf s}_{\text{tr}}=\left(\frac{{\bf s}}{c}\right)^{\frac{1}{\gamma}},$ (1)

where ${\bf s}$ is a sky model, $c$ is a constant, and $\gamma$ is a variable
parameter for which we explore a range of values, spanning from 1 to 30. The
value of $c$ is set to $2.96\times 10^{-5}$, corresponding to the maximum
value of image intensity identified empirically within the training set. This
method helps us to enhance the contrast of the image and mitigate the issues
arising from low signal strength. To map the variables to the range of
$[-1,1]$, we apply the following formula:

$\displaystyle{\bf x}_{0}=\frac{{\bf s}_{\text{tr}}-0.5}{0.5},$ (2)

where ${\bf x}_{0}$ is considered as a target image for the DDPM system.

The transformation in Equation 1 is not linear, and so it does not preserve
distances, which means that small relative values in the transformed image
correspond to a huge gap between values in the corresponding original image.
We illustrate the influence of power $\gamma$ from Equation 1 in Figure 2. In
this visualization, we take one sky model of the generated set and we show how
the suggested normalization in Equation 1 changes the image as well as its
histogram. The choice of $\gamma=1$ corresponds to the case where only linear
transform is applied. As $\gamma$ grows, the sources become wider and
brighter, and the histogram becomes flatter, indicating a denser image.

In the experiments, dirty images $\bf z$ are obtained from aggregated UV grids
$\bf y$, where

$\displaystyle{\bf z}=\frac{\frac{{\mathcal{F}}^{-1}({\bf
y})}{c_{d}}-0.5}{0.5},$ (3)

where $\mathcal{F}^{-1}$ is a function that maps UV grids to dirty images. In
CASA, $\mathcal{F}^{-1}$ includes inverse Fourier and artefact correction. The
$0.5$ factor is used to map images approximately to the $[-1,1]$ range, as in
Equation 2.

Images after $\mathcal{F}^{-1}$ are already dense and contain values of the
order of $c_{d}=10^{-5}$ Jy. For these images, we apply a scaling operation by
dividing each image by a factor of $c_{d}$.

These preprocessing steps map the sky model and dirty images to be within
standard ranges of typical machine learning models. The described DDPM
estimates the sky model denoted $\hat{\bf x}_{0}$. In order to map variables
back to the original scale, we apply the inverse operations:

$\displaystyle\hat{{\bf s}}_{\text{tr}}=0.5\hat{\bf x}_{0}+0.5,$ (4)
$\displaystyle\hat{{\bf s}}=c\hat{{\bf s}}_{\text{tr}}^{\gamma},$ (5)

where $c$ is the same constant as in Equation 1.

###  3.2 Image reconstruction based on stochastic models

In this study, we address image reconstruction in radio astronomy, which is an
inherently ill-posed problem. Information is lost due to UV sampling,
gridding, and handling of noise coming from the instruments and atmosphere.
Given these limitations, the task of accurately estimating the initial sky
model from the measured data is naturally ambiguous. Mathematically, it means
that many estimates $\hat{{\bf s}}$ might correspond to the same measured data
$\bf z$.

To estimate the sky model, we must inject priors into the model. There are two
possible approaches to do this: deterministic methods and stochastic methods.
Deterministic methods —as in CLEAN— use hand-crafted priors, and traditional
machine learning models are based on learnable frameworks; both operate in a
straightforward way. They create a one-to-one mapping: they first find optimal
parameters based on the training data and then, once these parameters are
fixed, they produce a deterministic output for each unique input. In contrast,
stochastic methods allow one-to-many mapping, meaning they can generate
multiple probable solutions for a single input. This is achieved by
introducing random variables or model-initiation noise into the model during
the computation of the output. Therefore, the output depends on both the input
data and the realization of the model initiation noise. By varying this model-
initiation noise for fixed model parameters, we can get different outputs for
the same input data. All such outputs are the most probable to align with the
target data manifold and correspond to the given input data (observations).
This inherent ambiguity is useful, as by inspecting different reconstructions
we can quantify the uncertainty on the predictions.

Figure 3: Schematic representation of diffusion process. The forward process
consists in progressively adding noise to the sky model. The inverse process
starts with model noise, which is progressively refined with a conditioning on
the dirty image.

###  3.3 Denoising diffusion probabilistic models

Figure 4: Generalized pipeline of proposed method. The top branch represents
inference. The DDPM model acts as a stochastic image-to-image translator,
transforming the conditioning $\bf z$ and the model noise $\bm{\epsilon}$ into
the target image $\hat{\bf x}_{0}$ and corresponding $\hat{{\bf s}}$. The
predicted sky model is then processed by the localization algorithm,
outputting the corresponding coordinates and fluxes. The bottom branch
represents the simplified measurement process. Figure 5: Schematic
representation of training pipeline for time step $t$. The preprocessed sky
model ${\bf x}_{0}$ is corrupted by the DDPM model noise $\bm{\epsilon}$. The
U-Net, given the corrupted target image $\mathbf{x}_{t}$ and conditioned with
the dirty image $\bf z,$ tries to estimate the model noise
$\hat{\bm{\epsilon}}$ to denoise the $\mathbf{x}_{t}$. The model is trained to
minimize the distance between ${\bm{\epsilon}}$ and $\hat{\bm{\epsilon}}$.

We use the image-to-image DDPM Palette Saharia et al. (2022). This diffusion
operates in an iterative manner, as visualized in Figure 3. There are two
directions: forward and inverse. During the forward direction stage, we add
model noise $\bm{\epsilon}$ to the data. This noise, an internal state
parameter of the model, is distinct from the astronomical noise arising from
antennas, the atmosphere, and so on. In DDPMs, this model noise is generally
sampled from a zero-mean Gaussian distribution. During the inverse direction
stage, the aim is to remove the added noise using a trained denoiser. This
denoiser is a crucial component of the denoising diffusion process and, in
this case, is represented by a U-Net (Ronneberger et al., 2015) with self-
attention layers (Vaswani et al., 2017). The parameters of U-Net are trained
to predict the model noise added to the target image. The iterative process of
applying this U-Net multiple times for progressive denoising is the function
of the DDPM.

At the forward process stage for each iteration $t\in\\{1,T\\},$ the added
model noise has its own variance, which we denote ${\alpha}_{1}$ to
${\alpha}_{T}$, respectively, where each ${\alpha}_{t}\forall t$ is between 0
and 1, and each subsequent value is greater than the previous one
${\alpha}_{T}>{\alpha}_{T-1}>...>{\alpha}_{1}$. Variances are chosen in such a
way that, after $T$ steps, we have a meaningless image corresponding to the
noise, essentially a sample from the Gaussian distribution. In this case, the
reverse process can be initiated from the noise.

Mathematically, the forward process can be written as

$\mathbf{x}_{t}=\alpha_{t}\mathbf{x}_{t-1}+\sqrt{1-{\alpha}_{t}}\bm{\epsilon},$
(6)

where $\mathbf{x}_{t}$ is an image produced by DDPM at the time step $t$, and
$\bm{\epsilon}$ is the added model state noise.

Ho et al. (2020) show that $\mathbf{x}_{t}$ can be rewritten as

$\mathbf{x}_{t}=\sqrt{\bar{\alpha}_{t}}\mathbf{x}_{0}+\sqrt{1-\bar{\alpha}_{t}}\bm{\epsilon},$
(7)

where ${\bf x}_{0}$ is a target data sample and
$\bar{\alpha}_{t}:=\prod_{s=1}^{t}\alpha_{s}$.

At the inverse process stage, we train the model to remove the added internal
noise (model noise). After $T$ steps, we assume that the model will encode a
target data manifold to itself, making it possible to generate data samples in
$T$ steps from a randomly initiated model state $\bm{\epsilon}$. In our
implementation, each step is conditioned on $\bf z$ by concatenating it to the
input. $\bf z$ represents the dirty image obtained from the aggregated UV grid
${\mathbf{y},}$ as shown in Figure 4 (bottom branch).

The training stage is visualized in Figure 5. The subset of 5 082 images is
used. During training, we uniformly select time steps $t\in[0,1000]$, which
define the noise variance for each instance. The target image (sky model)
$\mathbf{x}_{0}$ is first normalized, and then we generate an intermediate
$\mathbf{x}_{t}$ using Equation 7 and the model noise $\bm{\epsilon}$. This
$\mathbf{x}_{t}$ is concatenated with the normalized dirty image $\mathbf{z}$,
creating a two-channel tensor fed into the U-Net, which serves as our
denoiser. The U-Net architecture, which is adapted from Saharia et al. (2022),
has four stages of downsampling and four stages of upsampling. The time step
is encoded using two linear layers with SiLU activation functions. The
denoiser consists of residual blocks, which includes Group Normalization (with
32 groups and an epsilon value of $10^{-5}$), and Convolutional Layers with
kernel size 3. The model uses self-attention mechanisms. The U-Net accepts a
two-channel input and produces a one-channel output. This model contains
approximately 95M parameters. The primary training objective is to minimize
the $\ell_{2}$ distance between the actual model noise $\bm{\epsilon}$ and the
noise estimated by the U-Net
$\hat{\bm{\epsilon}}=g_{\bm{\theta}}(\mathbf{x}_{t},\mathbf{z})$, as defined
in Equation 8.

$\displaystyle{\mathcal{L}^{DDPM}(\bm{\theta})=\mathbb{E}_{t,{\bf z},{\bf
x}_{0},\bm{\epsilon}}\left[\left\|\bm{\epsilon}-g_{\bm{\theta}}\left(\sqrt{\bar{\alpha}_{t}}{\bf
x}_{0}+\sqrt{1-\bar{\alpha}_{t}}\bm{\epsilon},\bf
z,t\right)\right\|_{2}^{2}\right],}$ (8)

where ${\bf x}_{0}$ is a target image (sky model), $\bf z$ is a conditioning
(dirty image), $\bm{\epsilon}$ is Gaussian zero-mean unit variance noise
$\mathcal{N}({\mathbf{0}},{\mathbf{I}})$, $g_{\bm{\theta}}$ is conditional
DDPM, and $\bar{\alpha}_{t}=\prod_{s=1}^{t}\alpha_{s}$ is the noise variance
parameter, all at step $t,$ as in Equation 7. We train this U-Net over 237
epochs with a batch size of 24, and we use the Adam optimizer with a learning
rate of $10^{-4}$. We select the best model based on the loss on the
validation set consisting of 2541 images. The training process was completed
in approximately 16 hours on a NVIDIA RTX 3090 GPU.

The diffusion model can be seen as a stochastic image-to-image translator
whose goal consists in transforming dirty images $\bf z$ into a sky model
$\hat{\bf x}_{0}$. Therefore, dirty images $\bf z$ can also be referred to as
the input $\bf z$ to the DDPM alongside the model noise $\bm{\epsilon}$.

At the inference step, we begin with the Gaussian noise (noise of the model
$\bm{\epsilon}$) and the dirty image $\bf z$. This stage is not merely a
consecutive application of the denoiser to $\bm{\epsilon}$ but rather a
sophisticated iterative process that progressively refines the estimate of the
original sky model, $\hat{\bf x}_{0}$, starting with $\bm{\epsilon}$ through
$\mathbf{x}_{t}$. Each step of the DDPM involves predicting and then
subtracting the added noise from the current image $\mathbf{x}_{t}$ each time
concatenated with dirty image conditioning $\bf z$. After removing the
predicted noise from $\mathbf{x}_{t}$ we get an estimate of $\mathbf{x}_{0}$.
Then, a new, smaller-variance Gaussian noise is introduced to get
$\mathbf{x}_{t-1}$, which is then estimated and removed in the subsequent
step. Due to computational limitations, this process is repeated for a total
of 250 steps instead of 1000; it gradually reduces the noise and converges
toward an accurate reconstruction of $\mathbf{x}_{0}$. Subsequently, with
Equations 4 and 5, $\hat{\bf x}_{0}$ is mapped back to the original range of
values $\hat{{\bf s}}$. Finally, we apply the extraction algorithm to get the
properties of the sources from the predicted sky model. The process is shown
schematically in Figure 4 (top branch).

###  3.4 Source localization

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Photutils algorithm behavior on sky models for closely situated
sources.Behavior of the Photutils localization algorithm on sky models. Panels
(a) and (d) show the sky model; (b) and (e) depict images annotated with the
corresponding catalog sources; and (c) and (f) present images annotated with
sources detected using the chosen photutils algorithm. As demonstrated in
these images, the Photutils algorithm occasionally fails to deblend closely
situated sources. This happens for less than 1% of the sources present in the
dataset.

The locations, shapes, and fluxes of the sources can be restored from the
predicted sky model. While PyBDSF is commonly used for this task, it can only
be applied to images containing background noise. As the target sky models do
not contain any noise—just point sources—we find this extra step of adding
noise to execute PyBDSF destructive. This is especially true considering the
relatively straightforward nature of the sky models, which consist solely of
Gaussian sources with zero background. These challenges prompted us to explore
an alternative method of source localization.

We use an algorithm based on watershedding. We take an implementation of this
algorithm from the Photutils library by Bradley et al. (2016) to extract the
properties of point sources from predicted sky models. To validate the usage
of this algorithm for source localization, we evaluate its performance on all
sky models of the simulated dataset. The achieved scores are close to perfect,
as shown in Section 4. The algorithm fails when sources are very close to each
other (see Figure 6). Given the small amount of data for which it fails, we
use it as a final post-processing step of the pipeline.

The localization algorithm can be deployed directly on a single output from
the diffusion model; however, we can also use the stochasticity and run the
model with different realizations of model state noise to get multiple sky
model predictions. This can potentially enhance the model’s performance, but a
challenge arises in aggregating these multiple predictions. In this study, we
investigate two aggregation techniques to derive a final image from multiple
predictions: the mean and the median. These methods are applied at the pixel
level, where we assess the value of each pixel across all predictions and
compute either its mean or median. This process is repeated for every pixel,
resulting in the final output image.

####  3.4.1 Aggregation techniques

As mentioned above, DDPM can output multiple predictions for one phase of
conditioning. We investigate two ways to aggregate the produced predictions
via aggregate-detect and detect-aggregate strategies. The aggregate-detect
technique means that first we aggregate the images and then we localize the
sources. The detect-aggregate method means that we first localize sources from
each prediction and then we aggregate the sources into the final table.

For the aggregate-detect strategy, we study the mean and median as an
aggregation:

$\text{{Mean}}(h,w)=\frac{1}{n}\sum_{k=1}^{n}{\hat{\bf x}_{0}}^{k}(h,w),$ (9)

$\text{{Median}}(h,w)=\text{{median}}\left({\hat{\bf
x}_{0}}^{1}(h,w),{\hat{\bf x}_{0}}^{2}(h,w),\ldots,{\hat{\bf
x}_{0}}^{n}(h,w)\right),$ (10)

where $n$ is the total number of predicted images for a given $\bf z$ (20 in
this case), ${\hat{\bf x}_{0}}^{k}$ represents the $k$-th predicted image, and
$(h,w)$ represent the pixel coordinates.

The detect-aggregate method consists in aggregating the predicted sources from
multiple runs with the same conditioning. In this scenario, the source-
extraction algorithm is executed on each image independently. Subsequently,
the predictions are merged together; two predictions are considered to concern
the same source if their coordinates (right ascension (RA) and Declination
(Dec)) are within a radius of $r=5\times 10^{-5}$ degrees. In the final
aggregated output, each source is listed with its average coordinates and
characteristics (such as flux), along with the corresponding standard
deviation computed across the multiple predictions.

Additionally, we calculate the number of images in which each source is
detected, a value we term the “reliability” of the source. For instance, if a
source is detected in all 20 images out of 20, it is assigned a reliability of
100%; but if it is only detected in 10 images, the associated reliability is
50%. A reliability threshold can be set to determine which sources should be
reported. Intuitively, ambiguous sources will not appear in all predictions,
because they can also be classified as background noise. This should be
particularly true for very faint sources.

We show in Section 4 that as soon as the reliability threshold grows,
completeness falls, which means that the model starts to miss an increasing
number of sources (while purity grows). Reliability is expected to be a proxy
for the sensitivity of the proposed algorithm to faint sources.

###  3.5 Evaluation metrics

The task of localization in the described pipeline is divided into two stages:
sky model reconstruction from dirty images and source localization. We
evaluate the proposed method on both.

####  3.5.1 Reconstruction metrics

To evaluate how well the trained DDPM predicts the sky model from dirty image
$\bf z,$ we use metrics from natural image analysis such as peak S/N, SSIM,
and MSE, which are described below.

The mean squared error (MSE) quantifies the average squared difference between
the predicted and original pixel values:

$\text{MSE}=\frac{1}{H\times W}\sum_{h=1}^{H}\sum_{w=1}^{W}\biggl{(}\hat{\bf
x}_{0}(h,w)-{\bf x}_{0}(h,w)\biggr{)}^{2},$ (11)

where $H\times W$ represents the image size, ${\bf x}_{0}$ corresponds to the
target image, and $\hat{\bf x}_{0}$ represents the predicted image. A lower
MSE indicates better performance.

The peak signal-to-noise ratio (peak S/N) measures the quality of the
predicted image by comparing it to the target image; it is a calculation of
the ratio of the maximum possible pixel value MAX to MSE between the predicted
and original images:

$\text{peak S/N}=20\cdot\log_{10}(\text{MAX})-10\cdot\log_{10}(\text{MSE}).$
(12)

A higher peak S/N indicates better performance.

The structural similarity index (SSIM) introduced by Wang et al. (2004)
assesses the similarity between the predicted $\hat{\bf x}_{0}$ and target
image ${\bf x}_{0}$ based on their structural information:

${\mathrm{SSIM}}({\bf x}_{0},\hat{\bf x}_{0})=\frac{(2\mu_{{\bf
x}_{0}}\mu_{\hat{\bf x}_{0}}+c_{1})(2\sigma_{{\bf x}_{0}\hat{\bf
x}_{0}}+c_{2})}{(\mu_{{\bf x}_{0}}^{2}+\mu_{\hat{\bf
x}_{0}}^{2}+c_{1})(\sigma_{{\bf x}_{0}}^{2}+\sigma_{\hat{\bf
x}_{0}}^{2}+c_{2})},$ (13)

where $\mu_{{\bf x}_{0}}$ is the global mean value of ${\bf x}_{0}$ of size
$H\times W$; $\mu_{\hat{\bf x}_{0}}$ is the global mean value of $\hat{\bf
x}_{0}$ of size $H\times W$; $\sigma_{{\bf x}_{0}}$ is the global variance of
${\bf x}_{0}$; $\sigma_{\hat{\bf x}_{0}}$ is the global variance of $\hat{\bf
x}_{0}$; $\sigma_{{\bf x}_{0}\hat{\bf x}_{0}}$ is the covariance between ${\bf
x}_{0}$ and $\hat{\bf x}_{0}$; and $c_{1}=10^{-4}$ and $c_{2}=9\times 10^{-4}$
are two constants chosen to stabilize the division depending on the dynamic
range. The SSIM value ranges between 0 and 1, where a higher value indicates
better performance. The SSIM implementation used here is used from the skimage
package by Van der Walt et al. (2014).

####  3.5.2 Localization metrics

To evaluate the performance of the source localization task, we use purity and
completeness metrics following the same protocol as in Taran et al. (2023).
These metrics are essentially the same as precision and recall but we define
true positives (TPs), false positives (FPs), and false negatives (FNs)
differently. The predicted source position is considered a TP if its (RA, Dec)
coordinates are within a radius of $r=5\times 10^{-5}$ degrees from the true
source position, which is known from the simulation. A predicted source is a
FP if there is no true source within the same radius $r$. A FN is observed if
there are no predicted sources within the radius around a true source.

Purity, or precision, measures the accuracy of the predicted sources by
computing the ratio of the number of detected true sources (TP) to the total
number of predicted sources. It can be calculated using the formula:

$\text{Purity}=\frac{TP}{TP+FP}.$ (14)

Completeness, or recall, quantifies the completeness of the predictions by
calculating the ratio of the number of detected true sources (TP) to the total
number of true sources present in the data:

$\text{Completeness}=\frac{TP}{TP+FN}.$ (15)

In addition to purity and completeness, we also compute the F1-score, which is
their harmonic mean:

$\text{F1-score}=2\cdot\frac{\text{Purity}\cdot\text{Completeness}}{\text{Purity}+\text{Completeness}}.$
(16)

The F1-score ranges between 0 and 1, where a higher value indicates better
performance in terms of both purity and completeness.

####  3.5.3 Source characterization metrics

In the context of source characterization tasks, the goal is to predict the
flux values, which represent the integral brightness of the radio source. The
accuracy of the flux estimations is measured by calculating the proportion of
points for which the flux has been accurately predicted. We refer to a flux
prediction as accurate if the predicted value varies from the actual value
(established through the simulations) by a quantity equal to or less than the
noise amplitude. In the experiments, the noise mean square is set to 50 mJy.
Consequently, we introduce a metric known as the ”fraction”, which is defined
as the ratio of the number of points where the flux was accurately predicted
to the total number of sources in the simulation.

##  4 Results

###  4.1 Image reconstruction

The first stage of the proposed pipeline is a reconstruction of a sky model
from a dirty image.

####  4.1.1 Diffusion predictions

Dirty image

Sky model

Median

Mean

Std

Outputs

Figure 7: Diffusion model workflow: from dirty image to models outputs and
uncertainty estimation. Top row: (1) Original dirty image (conditioning), with
true sources circled in white. (2) The sky model showing sources detected by
the Photutils localization algorithm. (3) The post-processed predicted sky
model, aggregated using the median from 20 network outputs, with Photutils-
localized sources. (4) Another version of the post-processed predicted sky
model, this time aggregated pixel-wise using the mean from network outputs,
again with Photutils-localized sources. (5) The standard deviation of multiple
predictions before the post-processing step. We refer to this latter as the
uncertainty estimation; it is unitless, with brighter colors (like yellow)
indicating greater uncertainty. Rows 2-5 present 20 distinct outputs of the
diffusion model for the corresponding dirty image, each with a different model
noise initialization.

The model operates in a stochastic manner, producing multiple outputs for the
same conditioning image. This stochasticity provides advantages when dealing
with sources that exhibit low S/N. By running the model multiple times (20
runs in the experiments), we expect the faint sources to appear in the
fraction of the predictions in the same positions. It is important to note
that each run takes approximately 80 seconds per image on GPU RTX 3090. The
execution time can be further optimized using different generative techniques,
which will be the subject of future investigations.

Figure 7 shows an example of the model output. In the top row, we present
samples from the generated data set: the first image is the conditioning
(dirty image) and the second one is the corresponding sky model. The next
three outputs represent the aggregated mean, the aggregated median, and the
standard deviation of the predictions. The raw predictions of the model are in
the rows 2-5. We note that there exists a particularly faint source that the
proposed model fails to detect in aggregated images. However, in some single
runs the source is present. This leads to huge values of standard deviation in
this area, which shows that the model is uncertain about its predictions and
that further investigation is needed.

This decision to select 20 runs per dirty image is motivated by the observed
behavior of reconstruction metrics in Figure 8 for the median aggregation
process across multiple images. We notice an improvement in the metrics with
an increase in the number of aggregated predicted images. However, this
improvement comes at the expense of additional computation time. At the point
of $n=20$, we assume that further increases do not substantially improve the
metrics. Consequently, we identify $n=20$ as an effective compromise between
improved performance and computational efficiency.

Figure 8: Reconstruction metrics with respect to number of runs $n$ of
diffusion model for median aggregation. Table 1: Reconstruction results: MSE,
peak S/N, and SSIM metrics computed on preprocessed sky models for different
diffusion settings.

Normalization | Aggregation | Metrics ($\pm$ uncertainty)
---|---|---
MSE | peak S/N | SSIM
($\times 10^{-5}$) |  |
$\gamma=1$ | single run | $5.276\pm 0.138$ | $38.031\pm 1.152$ | $0.977\pm 0.002$
mean | $4.232\pm 0.333$ | $\mathit{42.243\pm 9.475}$ | $\mathit{0.998\pm 0.001}$
median | $\mathit{3.941\pm 0.354}$ | $41.587\pm 9.156$ | $\mathit{0.998\pm 0.001}$
$\gamma=2$ | single run | $1.492\pm 0.146$ | $38.999\pm 0.097$ | $0.830\pm 0.003$
mean | $\mathit{1.189\pm 0.604}$ | $41.251\pm 0.206$ | $0.853\pm 0.003$
median | $1.190\pm 0.605$ | $\mathit{42.148\pm 0.260}$ | $\mathit{0.882\pm 0.003}$
$\gamma=10$ | single run | $1.668\pm 0.171$ | $36.784\pm 0.464$ | $0.814\pm 0.004$
mean | $\mathit{1.183\pm 0.592}$ | $38.008\pm 2.263$ | $0.835\pm 0.005$
median | $\mathit{1.185\pm 0.586}$ | $\mathit{38.386\pm 2.550}$ | $\mathit{0.862\pm 0.005}$
$\gamma=20$ | single run | $2.264\pm 0.242$ | $35.105\pm 0.793$ | $0.892\pm 0.006$
mean | $\mathit{1.431\pm 0.645}$ | $37.457\pm 4.474$ | $0.929\pm 0.011$
median | $1.442\pm 0.643$ | $\mathit{37.651\pm 4.807}$ | $\mathit{0.955\pm 0.008}$
$\gamma=30$ | single run | $2.820\pm 0.313$ | $32.944\pm 1.086$ | $0.985\pm 0.002$
mean | $\mathit{1.553\pm 0.750}$ | $\mathit{34.844\pm 4.585}$ | $0.985\pm 0.008$
median | $1.569\pm 0.735$ | $34.693\pm 4.906$ | $\mathit{0.988\pm 0.007}$

Note: We varied the power root used in the normalization during training
(Equation 1) and the aggregation method (mean, median over 20 runs, or single
run). Single run signifies that the DDPM is run once for a given condition to
get the predicted sky model. We evaluate performance on each run individually
and then calculate the average score over all 20 runs. The error is calculated
by taking the standard deviation over all outputs for the test set for
aggregation methods, and over all runs to get a single run score.

The reconstruction metrics are presented in Table 1. These metrics are
calculated on the direct outputs of the diffusion model; that is, the
normalized sky model. We experiment with different values for the power root
$\gamma$ in the normalization function applied to the sky model (Equation 1).

As the original range of values within the images is significantly low,
metrics calculated within it would be exceedingly high, rendering them
uninformative. We therefore use normalized values. However, as discussed
above, higher root powers lead to lower contrast in the images. In this case,
the contrast is negatively correlated with the reconstruction metrics –-less
contrast in images tends to result in better reconstruction metrics for the
same model performance. This relationship is evident in Table 1, where lower
powers yield superior results compared to higher powers.

Aggregating multiple outputs from various runs outperforms the output from a
single run within the same normalization. This demonstrates that the inherent
stochasticity of the process assists in detecting fainter sources, which
indicates the value of using statistics across multiple runs to enhance
reconstruction. Through experimentation, we find that the mean and median
aggregation methods deliver very similar results.

An essential takeaway from this analysis is that reconstruction metrics may
not be the most adequate indicators for assessing the quality of astronomical
data. In subsection 4.2, we explore the localization metrics and discuss their
performance for the same settings.

####  4.1.2 Uncertainty estimation in the image space

To quantify the uncertainty associated with the sky model predictions, we use
the concept of uncertainty estimation in the image space. We define
$\mathbf{x}_{\text{uncertainty}}$ using the following formula:

$\mathbf{x}_{\text{uncertainty}}(h,w)=\sqrt{\frac{1}{n}\sum_{j=1}^{n}\left(\hat{\bf
x}_{0}^{j}(h,w)-\frac{1}{n}\sum_{k=1}^{n}\hat{\bf
x}_{0}^{k}(h,w)\right)^{2}},$ (17)

where $n$ represents the number of runs ($n=20$), $\hat{\bf x}_{0}^{j}$
denotes a prediction of the sky model for the $j$-th run, and $(h,w)$ are the
pixel coordinates. This formula denotes the standard deviation visualized in
the top right of Figure 7. This latter displays the variability in the
predictions across multiple runs. Intuitively, if the model consistently
identifies the same sources, the uncertainty is low and we can be certain
about their presence. Conversely, if the sources appear and disappear from the
predictions across runs, the variance in the image space will be high,
indicating a region of the sky that requires further investigation. The
variance in the predictions can be juxtaposed with the intensity values
estimated in the mean-aggregated image to determine the degree of uncertainty
associated with a particular region. This can inform the final decision about
whether a given region exhibits significant uncertainty, thereby guiding
further analysis.

###  4.2 Source localization

After obtaining a predicted sky model, we try to identify the sources present.
We evaluate the performance of our source localization using three metrics:
purity, completeness, and F1-score, as detailed in Table 2.

The F1-score is a particularly effective measure of model performance because
it simultaneously considers both purity (the model’s ability to limit FPs) and
completeness (the model’s ability to identify all real sources in the
dataset). We therefore consider F1 as the main indicator for choosing the best
model.

We provide the scores for a Photutils method of source localization applied
directly to the simulated sky model given by CASA (last line in Table 2). This
reference reveals the level of error that comes from the chosen parameters of
the algorithm from Photutils.222The Photutils parameters are: kernel_x_size=3,
kernel_y_size=3, sigma_clipped_stats_sigma=2. detect_npixels=10,
deblend_sources_npixels=10, std_const=120.  These parameters stay the same
during all experiments, which ensures that any differences in outcomes can be
attributed to the models themselves, rather than to variations in the
processing steps.

We draw a comparison with the model presented by Taran et al. (2023). The
distinguishing feature of their model is an additional compression step, which
retains only 1450 fixed positions in the UV plane and operates directly in the
frequency domain. Contrasting this with the proposed DDPM, which employs 250
passes, the model in Taran et al. (2023) accomplishes its task in a single
pass. This characteristic exemplifies the common trade-off between accuracy
and computational efficiency. The model in Taran et al. (2023) performs better
than PyBDSF (Mohan and Rafferty, 2015), which currently thought to be state of
the art. Consequently, we chose to compare our results to those of obtained
using the model in Taran et al. (2023). As can be seen from Table 2, all the
models achieve very high F1-scores exceeding 95 %, compared to 81 % achieved
by the Taran et al. (2023) model.

For the images with lower root functions $\gamma\in[1,2]$, the best
aggregation methods are aggregate-detect ones. The outputs for the model
trained with such normalizations are noisy, which is due to the sparsity of
the images. Aggregation on the image level averages the noise in each
prediction, producing a higher-quality output. However, as $\gamma$ gets
larger and the prediction quality rises, the winning method becomes detect-
aggregate.

The optimal model employs $\gamma=2$. We illustrate the highest localization
scores relative to various $\gamma$ values in Figure 9. The $\gamma=2$ case
demonstrates the effectiveness of the proposed normalization technique
compared to the $\gamma=1$ case. However, as $\gamma$ increases, the
performance falls because of excessive sharpening of the blob borders, as seen
in Figure 2.

In Figure 10, we plot the scores for $\gamma=2$ against the number of runs. As
$n$ increases, the F1-score saturates, reaching a plateau at $n=5$ with no
significant further improvement. This behavior contrasts with the
reconstruction scenario, where the metrics continue to improve with increasing
$n$. This occurs for reconstruction metrics because the more images we
generate, the more errors are corrected when averaging, leading to smooth
improvement. However, localization metrics operate differently; their primary
focus is on the position of the flux. Consequently, fewer images are required
for the flux to appear in the final prediction, and adding more images does
not significantly modify the chosen localization metrics. Essentially, these
metrics are only concerned with whether the source was detected or not, which
is why they exhibit a different behavior.

Figure 9: Best scores for purity, completeness, and F1 among all aggregations
with respect to different root power $\gamma$ used in normalization during
training. Figure 10: Localization metrics with respect to number of runs $n$
of diffusion model.

###  4.3 Uncertainty in localization

The detect-aggregate technique confers an advantage, in that it suggests a
clear way to find uncertainty on the estimations. We compute the percentage of
images where the source is present and then we aggregate those source
characteristics through averaging. The fraction of output images where the
source is present is interpreted as a measure of source reliability: in this
study, the source is not retained if the reliability is less than 30%.

To illustrate the rationale behind the choice of the 30% threshold, we plot
the F1-score, purity, and completeness against the reliability threshold for
retaining sources in Figure 11. All curves except the one with root power
$\gamma=1$ reach saturation around 30%, have a flat section after, and start
to diminish in the end. Purity grows as we demand that the source be more
consistent in the images, and completeness falls as we eliminate an increasing
number of sources that appear in an insufficient number of images.

We plot the dependence of reliability on the S/N of the sources in Fig. 12. We
find that most of the sources with low reliability have low-S/N values. This
demonstrates that the stochasticity of the DDPM model mostly affects the
low-S/N sources, and they therefore have lower reliability. Consequently, for
sources detected with low S/N using the proposed model, reliability can be
used as an additional metric to select more robust sources. Choosing a
reliability threshold of 30% is a well-balanced trade-off for the simulated
data, allowing the model to have few FPs and FNs even at relatively low S/N.
This is discussed in further detail in Section 4.4.

Figure 11: Dependence of F1, purity, and completeness scores on the
reliability parameter for all models. Reliability is the fraction of the 20
images in which a source is present. The closest to optimal value we find is
30%, which means each source should be present at least in 6 images out of 20.
Table 2: Purity, completeness, and F1-score.

Model | Normalization | Input | Aggregation | Metrics
---|---|---|---|---
Purity | Completeness | F1
Diffusion | $\gamma=1$ | dirty image | single run | 100.00 | 44.97 | 62.04
detect-aggregate | 100.00 | 53.52 | 69.72
mean | 99.60 | 93.11 | 96.24
median | 99.83 | 90.84 | 95.12
Diffusion | $\gamma=2$ | dirty image | single run | 97.45 | 96.71 | 97.08
detect-aggregate | 98.33 | 97.88 | 98.10
mean | 97.88 | 97.80 | 97.84
median | 99.30 | 97.05 | 98.16
Diffusion | $\gamma=10$ | dirty image | single run | 98.66 | 94.66 | 96.62
detect-aggregate | 98.90 | 95.69 | 97.27
mean | 98.51 | 95.23 | 96.84
median | 99.29 | 94.78 | 96.98
Diffusion | $\gamma=20$ | dirty image | single run | 98.26 | 94.43 | 95.81
detect-aggregate | 99.52 | 94.70 | 97.05
mean | 99.52 | 93.64 | 96.49
median | 99.52 | 93.79 | 96.57
Diffusion | $\gamma=30$ | dirty image | single run | 98.58 | 93.46 | 94.93
detect-aggregate | 99.20 | 94.10 | 96.58
mean | 98.80 | 93.41 | 96.03
median | 98.81 | 94.25 | 96.47
PyBDSF | - | clean image | - | 72.18 | 20.82 | 32.31
Taran et al. (2023) | - | reduced uv-samples | - | 91.02 | 74.14 | 81.72
Photutils localization | - | sky model | - | 99.70 | 99.10 | 99.40

Note: The localization is performed on predictions that have been renormalized
to be within the same range as the original sky model. The model
configurations vary in terms of normalization power used during training, and
also in terms of the aggregation method employed. For the normalization with
high root powers $\gamma,$ the optimal aggregation method is detect-aggregate.
This approach offers a good trade-off between purity and completeness,
resulting in the best F1-score. For the lower $\gamma,$ the best approach
turned out to be aggregate-detect. We compare the results with Mohan and
Rafferty (2015) (PyBDSF) and Taran et al. (2023). We run the Photutils
algorithm directly on sky models to quantify the error coming from this step,
denoted as Photutils localization.

###  4.4 Performance vs S/N

We study the performance of the model as a function of S/N. We use a
definition of normalized injected flux from Béthermin et al. (2020):

$f_{\text{norm}}=\frac{\text{total
flux}}{\sigma_{\text{noise}}}\frac{b_{\text{min}}\cdot
b_{\text{maj}}}{{\sqrt{b^{2}_{\text{min}}+s^{2}_{\text{min}}}}\cdot\sqrt{b^{2}_{\text{maj}}+s^{2}_{\text{maj}}}},$

where $b_{\text{min}}$ and $b_{\text{maj}}$ are the minor and major axes of
the beam (0.82” and 0.89” respectively) and $s_{\text{min}}$ and
$s_{\text{maj}}$ are the minor and major axes of the source. The results are
shown in Figure 13 for the best aggregation method among the tested diffusion
models (see Table 2; the best are highlighted in bold for each model and
normalization). We conduct a comprehensive comparison with the model from
Taran et al. (2023). We also compare the results to those obtained by
Béthermin et al. (2020).

All of the models showcase a robust and significant performance boost,
particularly for low S/N, in comparison to previous methods. This improvement
is highlighted by the fact that all models reach a completeness score
exceeding 90% even at a low S/N of 2, which can be compared to the score of
50% found by Taran et al. (2023) for the same S/N. For S/N values higher than
3, all models reach close to 100% completeness, which means that they do not
miss any real sources.

The reliability assigned through the detect-aggregate method is anticipated to
depend on the S/N of the source, as demonstrated in Figure 12. We observe a
trend where sources with lower reliability tend to have lower S/Ns. However,
it is notable that the model is able to consistently reconstruct many sources,
even those with low S/N.

Figure 12: Signal-to-noise ratio versus reliability. Lower reliability values
are predominantly associated with sources that show lower S/Ns. However, the
model often localizes sources with low S/Ns with high confidence. Figure 13:
Completeness comparison with state-of-the-art methods. These state-of-the-art
methods outperform traditional algorithms such as PyBDSF, which is why the
comparison is limited to them. All diffusion models demonstrate significant
improvement at low S/N.

###  4.5 Flux estimation

Figure 14: Relative error in flux estimation compared to true flux for studied
aggregation techniques across images. The gray area indicates the 1$\sigma$
range. Per-pixel aggregation methods, such as the mean and median, yield
reliable estimates. Visually, the median provides the most appealing results
with the fewest outliers. All results presented here are for the model
employing $\gamma=2$. The outliers are discussed in Appendix A.

As the proposed model is trained to recover the full sky, we also attempted to
estimate the flux of the sources in addition to their locations. Here too, we
use the Photutils algorithm to estimate the flux of the sources in the
predicted sky models. The model predicts a normalized sky model. The
normalization is easily invertible, allowing us to estimate the fluxes
directly after renormalization of the model outputs. We present the results
for all aggregation methods in Figures 14 and 15. We plot the relative error
$\frac{\text{flux}_{\text{estimated}}-\text{flux}_{\text{true}}}{\text{flux}_{\text{true}}}$
on the y-axis and the true flux ($\text{flux}_{\text{true}}$) on the x-axis,
where $\text{flux}_{\text{estimated}}$ is the estimated flux from the
predicted sky model and $\text{flux}_{\text{true}}$ is known from the
simulation.

The gray area represents the allowable error range: the root-mean-squared
(rms) noise of the simulated data is $\sim$50 $\mu$Jy. If the flux estimation
error falls within this range, it indicates that the model’s estimations are
robust.

Figure 15: Relative error in flux estimation as compared to true flux,
aggregated across sources from each predicted image (detect-aggregate method).
The gray area indicates the 1$\sigma$ range. The error bars represent the
standard deviation of the flux estimates obtained from each prediction. The
color scale represents the reliability. The results showcased here are based
on the model that uses $\gamma=2$. The outliers are discussed in the Appendix
A.

In Figure 14, we compare different aggregation methods and examine their
impact on flux estimation. From this plot, it is clear that aggregation
substantially improves the estimates as it averages out errors across all
predictions, further highlighting the promise of using stochastic models.
Notably, the single run does not yield good results, indicating the importance
of per-pixel aggregation. As expected, median aggregation has fewer outliers
than mean aggregation. Aggregating through the mean is more efficient than
through detect-aggregate, making it a viable alternative.

The detect-aggregate method is shown in Figure 15. Here for the flux of each
source, we can also estimate the uncertainty in the flux by measuring the
standard deviation of the flux across several runs. We show this as error bars
on the fluxes in Figure 15. The reliability of each source is color-coded: in
the chosen color map, blue corresponds to sources detected in small fractions
of predictions, whereas red corresponds to sources that were consistently
detected across all of them, giving a reliability of 100 %. We notice that
there are very few outliers in the flux estimation. Interestingly, their
reliability scores are high. It is remarkable that the fluxes of the low-S/N
sources are estimated relatively accurately.

Table 3: Comparative results for the fraction of the detected sources whose
estimated flux is within 50 $\mu$Jy of the true flux.

Model | Normalization | Aggregation
---|---|---
Mean | Median | Single run | Detect-aggregate
CLEAN+PyBDSF | - | 0.57
Diffusion | $\gamma=1$ | 0.85 | 0.66 | 0.048 | 0.015
$\gamma=2$ | 0.96 | 0.97 | 0.88 | 0.97
$\gamma=10$ | 0.92 | 0.93 | 0.71 | 0.89
$\gamma=20$ | 0.82 | 0.83 | 0.53 | 0.72
$\gamma=30$ | 0.74 | 0.75 | 0.47 | 0.66

Note: The detect-aggregate method results are reported for a reliability
threshold of 30%. We compare our results with PyBDSF.

In Table 3, we compare the proposed method against CLEAN+PyBDSF, measuring the
fraction of sources with the correct flux estimation (whithin 50 $\mu$Jy)
among all sources detected. The method significantly outperforms PyBDSF
scores. We also notice that any aggregation improves the score compared to the
single run for $\gamma\geq 2$.

The model trained with $\gamma=1$ demonstrates low scores for flux estimation
in both single-run and detect-aggregate methods. This underperformance can be
attributed to the impulsiveness of the predicted images and the inability of
the network to effectively reconstruct them in a single run. In this case, the
detect-aggregate method also fails, as the errors from each individual image
are considerable. However, the aggregate-detect approach significantly
improves the scores, even for the case of $\gamma=1$. For both localization
and flux-estimation tasks, the best performance is achieved with the model
trained using the sky model normalization approach with $\gamma=2$ and median
aggregation method.

###  4.6 Sensitivity to the noise

The data set used for training the model is characterized by a noise level of
50 $\mu$Jy. During real observations, the noise level does not remain
constant. One of the sources of noise is the water vapor present in the
atmosphere. To ensure that the model is applicable in real settings, it is
vital that it performs well not only under the conditions it was trained on
but also when exposed to varying atmospheric conditions. To evaluate how the
performance of the model varies in response to different noise levels, we
conducted an analysis in which we alter the amount of water vapor. In order to
do this, we used the additional simulations with varying PWV parameter. We
tested the model on these data.

Figure 16: Impact of water vapor on reconstruction metrics. The MSE, SSIM, and
peak S/N values are normalized to highlight trends. As the amount of water
vapor in the atmosphere increases, the quality of the reconstruction
diminishes, evidenced by the decreasing SSIM and peak S/N, and the increasing
MSE.

The reconstruction results are shown in Figure 16. We see that all metrics
decrease as the noise increases. The model did experience different S/N during
training, which explains the relatively small decrease in the performance.

Figure 17: Impact of water vapor on localization metrics. The purity exhibits
a relative change of about 1%, while the completeness decreases significantly
by 4%. The model does not identify many FPs as the noise levels increase, but
its capability to accurately detect existing sources diminishes rapidly.
Notably, even with the highest noise levels, the completeness score
outperforms the previous scores on the test set with the low water vapor
noise. The dashed line is at 1.796 PWV units, which corresponds to the
training data. Figure 18: Influence of water vapor on fraction of point
sources whose predicted flux falls within noise amplitude from true flux. As
water vapor levels rise, this fraction significantly decreases. The dashed
line represents a PWV value of 1.796 units, which corresponds to the
conditions of the training data.

The localization results are illustrated in Figure 17, where we plot the
model’s completeness and purity as a function of the amount of water vapor
(and therefore noise level). From the figure, it is evident that as the amount
of water vapor increases, the model’s completeness decreases. This decline in
completeness indicates that the model is less sensitive to faint sources when
they are affected by noise.

In contrast, the purity of the model does not exhibit a significant decrease
as the amount of water vapor increases. The relative stability of purity
suggests that, even though the model fails to detect an increasing number of
true sources as the noise levels rise (number of FNs rises), it does not
produce many FPs by misinterpreting noise as actual sources.

However, changing the water vapor has a significant influence on flux
estimation. We present the results of these experiments in Figure 18. It is
evident that the best performance occurs at 1.796 PWV units, which corresponds
to the conditions of the training set. Among all the methods, detect-aggregate
proves to be the most robust to changes in noise. Nevertheless, all methods
show a significant decrease in performance as water vapor levels increase.
Training the model for specific water vapor values is needed to maintain
performance.

##  5 Discussion

The method introduced in this paper demonstrates state-of-the-art performance
in localizing point sources in “dirty” radio astronomy images, even under
low-S/Ns conditions. The results showcase the capability of the presented
model to outperform existing methods such as PyBDSF in both localization and
characterization of the sources, and particularly in estimating the flux of
the sources for the simulated dataset.

The primary limitation of the model is its reconstruction speed. Using an RTX
3090 GPU, the model executes 250 DDPM steps for each of the 20 realizations
and takes approximately one minute and eight seconds to produce a single
image. This process is conducted with a batch size of 20, allowing for the
parallel processing of 20 realizations simultaneously, thereby optimizing the
GPU’s capabilities. We note that running these 20 parallel realizations
requires a GPU with at least 8 GB of memory.This bottleneck presents an area
for future development, with a need to explore more efficient generation
techniques from diffusion models. Optimizing the speed of a DDPM is one of the
main problems for the engineering research community and several solutions
have already been suggested by Karras et al. (2022).

An additional direction of interest is to explore the transferability of the
proposed diffusion models. We would like to see whether or not a model trained
on data from one telescope can successfully be transferred to another.
Similarly, how well a DDPM model trained on simulated data can perform when
applied to real data remains unexplored. Also, the portability of the model to
SKA is unknown. These queries highlight potential limitations of the current
approach and we plan to address them in future work.

##  6 Conclusions

In this work, we introduce a new approach for interferometric image
reconstruction in astronomy using diffusion models, in particular the
conditional denoising diffusion probabilistic model (DDPM). The DDPM model is
trained, validated, and tested on simulated deep ALMA observations of extra-
galactic targets previously used by Taran et al. (2023). In the proposed
approach, we take the dirty images produced from simulated UV visibilities as
input with some preprocessing and return a set of different sky models. Given
the inherent stochasticity of DDPM models, the described model produces
multiple sky models for the same dirty image. We propose two different
approaches to aggregate the predictions: aggregate-detect and detect-
aggregate. The properties of the sources are extracted from predicted sky
models using the Photutils algorithm. We compare locations and flux values
from our DDPM model with the true sky model from the simulations. We devise a
simple diagnostic to measure the reliability of a source in predicted sky
models. We demonstrate that the DDPM model can be used to perform source
localization with high purity and completeness. Finally, we show that it can
also be used to accurately estimate source flux.

Our main results can be summarized as follows:

  * •

Aggregating the sky models from multiple runs of the described model using
mean or median techniques produces a good representation of the true sky
model. Producing a larger number of such runs ($\sim$ 20) produces aggregated
sky models with low MSE and high peak S/N and SSIM.

  * •

The proposed model further performs an even better image reconstruction if we
preprocess the true sky models while training with a power $\gamma=2$ in
Equation 1.

  * •

The DDPM-based approach provides a significant improvement in source
localization over traditional methods in terms of purity and completeness,
particularly at low S/N.

  * •

For S/N=2, the described model demonstrates 70% accuracy without any
normalization techniques, outperforming the previous state-of-art algorithm,
PyBDSF, which achieves 55% accuracy.

  * •

We find that $n=5$ runs is optimal, allowing the model to obtain an F1-score
of above 98 %.

  * •

We introduce a method to estimate the reliability of the predicted sources
using the inherent stochastic nature of DDPM.

  * •

We also propose a framework for estimating the fluxes of sources as the method
focuses on reconstructing the sky model. The performance surpasses that of
PyBDSF by 32% in terms of the fraction of sources for which the flux is
estimated within the noise level.

  * •

We also examine the influence of the normalization power, $\gamma$, of the
root transformation on the efficiency of the method. Our findings suggest that
training with $\gamma=2$ is optimal for both localization and flux estimation.

  * •

We study the impact of noise magnitude by varying water vapor levels. The
model’s completeness drops to 88%, yet its purity stays consistent, suggesting
that it identifies few false positives despite increased noise. Adapting the
model to specific water vapor values is crucial when using it for flux
estimation, though the proposed approach consistently outperforms PyBDSF.

In summary, we developed a novel DDPM-based method for source detection and
flux estimation that operates directly on dirty images. This approach offers
significant improvements over existing algorithms in source detection across
all S/Ns.

## Acknowledgements

M. Drozdova and V. Kinakh are supported by the _RODEM: Robust deep density
models for high-energy particle physics and solar flare analysis_ Sinergia
Project funded by the Swiss National Science Foundation, grant number
CRSII5-193716. O. Taran and O. Bait are supported by the _AstroSignals_
Sinergia Project funded by the Swiss National Science Foundation, grant number
CRSII5-193826. We also would like to thank the editor and referee for their
useful comments.

## References

  * Akiyama et al. [2019] Kazunori Akiyama, Antxon Alberdi, Walter Alef, Keiichi Asada, Rebecca Azulay, Anne-Kathrin Baczko, David Ball, Mislav Baloković, John Barrett, Dan Bintley, et al.  First m87 event horizon telescope results. ii. array and instrumentation.  _The Astrophysical Journal Letters_ , 875(1):L2, 2019.
  * Bean et al. [2022] Ben Bean, Sanjay Bhatnagar, Sandra Castro, Jennifer Donovan Meyer, Bjorn Emonts, Enrique Garcia, Robert Garwood, Kumar Golap, Justo Gonzalez Villalba, Pamela Harris, et al.  Casa, common astronomy software applications for radio astronomy.  _Publications of the Astronomical Society of the Pacific_ , 134(1041):114501, 2022.
  * Béthermin et al. [2020] M. Béthermin, Y. Fudamoto, M. Ginolfi, F. Loiacono, Y. Khusanova, P. L. Capak, P. Cassata, A. Faisst, O. Le Fèvre, D. Schaerer, J. D. Silverman, L. Yan, R. Amorin, S. Bardelli, M. Boquien, A. Cimatti, I. Davidzon, M. Dessauges-Zavadsky, S. Fujimoto, C. Gruppioni, N. P. Hathi, E. Ibar, G. C. Jones, A. M. Koekemoer, G. Lagache, B. C. Lemaux, C. Moreau, P. A. Oesch, F. Pozzi, D. A. Riechers, M. Talia, S. Toft, L. Vallini, D. Vergani, G. Zamorani, and E. Zucca.  The ALPINE-ALMA [CII] survey: Data processing, catalogs, and statistical source properties.  _Astronomy & Astrophysics_, 643:A2, oct 2020.  doi: 10.1051/0004-6361/202037649.  URL https://doi.org/10.1051%2F0004-6361%2F202037649.
  * Béthermin et al. [2020] M Béthermin, Y Fudamoto, M Ginolfi, F Loiacono, Y Khusanova, PL Capak, P Cassata, A Faisst, O Le Fèvre, D Schaerer, et al.  The alpine-alma [cii] survey: Data processing, catalogs, and statistical source properties.  _Astronomy & Astrophysics_, 643:A2, 2020.
  * Bishop [2007] Christopher M. Bishop.  _Pattern Recognition and Machine Learning (Information Science and Statistics)_.  Springer, 1 edition, 2007.  ISBN 0387310738.  URL http://www.amazon.com/Pattern-Recognition-Learning-Information-Statistics/dp/0387310738%3FSubscriptionId%3D13CT5CVB80YFWJEPWS02%26tag%3Dws%26linkCode%3Dxm2%26camp%3D2025%26creative%3D165953%26creativeASIN%3D0387310738.
  * Bouwens et al. [2022] R. J. Bouwens, R. Smit, S. Schouws, M. Stefanon, R. Bowler, R. Endsley, V. Gonzalez, H. Inami, D. Stark, P. Oesch, J. Hodge, M. Aravena, E. da Cunha, P. Dayal, I. de Looze, A. Ferrara, Y. Fudamoto, L. Graziani, C. Li, T. Nanayakkara, A. Pallottini, R. Schneider, L. Sommovigo, M. Topping, P. van der Werf, H. Algera, L. Barrufet, A. Hygate, I. Labbé, D. Riechers, and J. Witstok.  Reionization Era Bright Emission Line Survey: Selection and Characterization of Luminous Interstellar Medium Reservoirs in the z ¿ 6.5 Universe.  _ApJ_ , 931(2):160, June 2022.  doi: 10.3847/1538-4357/ac5a4a.
  * Bradley et al. [2016] Larry Bradley, Brigitta Sipocz, Thomas Robitaille, Erik Tollerud, Christoph Deil, Zè Vinícius, Kyle Barbary, Hans Moritz Günther, Azalee Bostroem, Michael Droettboom, et al.  Photutils: Photometry tools.  _Astrophysics Source Code Library_ , pages ascl–1609, 2016.
  * Clark [1980] B. G. Clark.  An efficient implementation of the algorithm ’CLEAN’.  _A &A_, 89(3):377, September 1980.
  * Connor et al. [2022] Liam Connor, Katherine L Bouman, Vikram Ravi, and Gregg Hallinan.  Deep radio-interferometric imaging with polish: Dsa-2000 and weak lensing.  _Monthly Notices of the Royal Astronomical Society_ , 514(2):2614–2626, 2022.
  * Conway et al. [1990] J. E. Conway, T. J. Cornwell, and P. N. Wilkinson.  Multi-frequency synthesis : a new technique in radio interferometrie imaging.  _MNRAS_ , 246:490, October 1990.
  * Cornwell and Evans [1985] T J Cornwell and KF Evans.  A simple maximum entropy deconvolution algorithm.  _Astronomy and Astrophysics (ISSN 0004-6361), vol. 143, no. 1, Feb. 1985, p. 77-83._ , 143:77–83, 1985.
  * Dhariwal and Nichol [2021] Prafulla Dhariwal and Alexander Nichol.  Diffusion models beat gans on image synthesis.  _Advances in neural information processing systems_ , 34:8780–8794, 2021.
  * Faisst et al. [2020] A. L. Faisst, D. Schaerer, B. C. Lemaux, P. A. Oesch, Y. Fudamoto, P. Cassata, M. Béthermin, P. L. Capak, O. Le Fèvre, J. D. Silverman, L. Yan, M. Ginolfi, A. M. Koekemoer, L. Morselli, R. Amorín, S. Bardelli, M. Boquien, G. Brammer, A. Cimatti, M. Dessauges-Zavadsky, S. Fujimoto, C. Gruppioni, N. P. Hathi, S. Hemmati, E. Ibar, G. C. Jones, Y. Khusanova, F. Loiacono, F. Pozzi, M. Talia, L. A. M. Tasca, D. A. Riechers, G. Rodighiero, M. Romano, N. Scoville, S. Toft, L. Vallini, D. Vergani, G. Zamorani, and E. Zucca.  The ALPINE-ALMA [C II] Survey: Multiwavelength Ancillary Data and Basic Physical Measurements.  _ApJS_ , 247(2):61, April 2020.  doi: 10.3847/1538-4365/ab7ccd.
  * Finke et al. [2021] Thorben Finke, Michael Krämer, Alessandro Morandini, Alexander Mück, and Ivan Oleksiyuk.  Autoencoders for unsupervised anomaly detection in high energy physics.  _Journal of High Energy Physics_ , 2021(6), jun 2021.  doi: 10.1007/jhep06(2021)161.  URL https://doi.org/10.1007%2Fjhep06%282021%29161.
  * Franco et al. [2018] M Franco, D Elbaz, M Béthermin, B Magnelli, C Schreiber, L Ciesla, M Dickinson, N Nagar, J Silverman, Emanuele Daddi, et al.  Goods-alma: 1.1 mm galaxy survey-i. source catalog and optically dark galaxies.  _Astronomy & Astrophysics_, 620:A152, 2018.
  * Gheller and Vazza [2022] Claudio Gheller and Franco Vazza.  Convolutional deep denoising autoencoders for radio astronomical images.  _Monthly Notices of the Royal Astronomical Society_ , 509(1):990–1009, 2022.
  * Gull and Daniell [1978] Stephen F. Gull and Geoffrey J. Daniell.  Image reconstruction from incomplete and noisy data.  _Nature_ , 272:686–690, 1978.  URL https://api.semanticscholar.org/CorpusID:4163124.
  * Ho et al. [2020] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  _Advances in neural information processing systems_ , 33:6840–6851, 2020.
  * Högbom [1974] J. A. Högbom.  Aperture Synthesis with a Non-Regular Distribution of Interferometer Baselines.  _A &AS_, 15:417, June 1974.
  * Högbom [1974] JA Högbom.  Aperture synthesis with a non-regular distribution of interferometer baselines.  _Astronomy and Astrophysics Supplement Series_ , 15:417, 1974.
  * Hotan et al. [2021] AW Hotan, JD Bunton, AP Chippendale, M Whiting, J Tuthill, VA Moss, D McConnell, SW Amy, MT Huynh, JR Allison, et al.  Australian square kilometre array pathfinder: I. system description.  _Publications of the Astronomical Society of Australia_ , 38:e009, 2021.
  * Jonas and MeerKAT Team [2016] J. Jonas and MeerKAT Team.  The MeerKAT Radio Telescope.  In _MeerKAT Science: On the Pathway to the SKA_ , page 1, January 2016.  doi: 10.22323/1.277.0001.
  * Karras et al. [2022] Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine.  Elucidating the design space of diffusion-based generative models.  _Advances in Neural Information Processing Systems_ , 35:26565–26577, 2022.
  * Le Fèvre et al. [2020] O Le Fèvre, M Béthermin, A Faisst, GC Jones, P Capak, P Cassata, JD Silverman, D Schaerer, L Yan, R Amorin, et al.  The alpine-alma [cii] survey-survey strategy, observations, and sample properties of 118 star-forming galaxies at 4¡ z¡ 6.  _Astronomy & Astrophysics_, 643:A1, 2020.
  * Liu et al. [2019] Daizhong Liu, Eva Schinnerer, Brent Groves, B Magnelli, Philipp Lang, S Leslie, E Jiménez-Andrade, DA Riechers, Gergo Popping, Georgios E Magdis, et al.  Automated mining of the alma archive in the cosmos field (a3cosmos). ii. cold molecular gas evolution out to redshift 6.  _The Astrophysical Journal_ , 887(2):235, 2019.
  * McMullin et al. [2007] Joseph P McMullin, BSDYWGK Waters, Darrell Schiebel, Wei Young, and Kumar Golap.  Casa architecture and applications.  In _Astronomical data analysis software and systems XVI_ , volume 376, page 127, 2007.
  * Mohan and Rafferty [2015] Niruj Mohan and David Rafferty.  Pybdsf: Python blob detection and source finder.  _Astrophysics Source Code Library_ , pages ascl–1502, 2015.
  * Narayan and Nityananda [1986] Ramesh Narayan and Rajaram Nityananda.  Maximum entropy image restoration in astronomy.  _Annual review of astronomy and astrophysics_ , 24(1):127–170, 1986.
  * Ramesh et al. [2022] Aditya Ramesh, Prafulla Dhariwal, Alex Nichol, Casey Chu, and Mark Chen.  Hierarchical text-conditional image generation with clip latents.  _arXiv preprint arXiv:2204.06125_ , 1(2):3, 2022.
  * Rau and Cornwell [2011] Urvashi Rau and Tim J Cornwell.  A multi-scale multi-frequency deconvolution algorithm for synthesis imaging in radio interferometry.  _Astronomy & Astrophysics_, 532:A71, 2011.
  * Ronneberger et al. [2015] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18_ , pages 234–241. Springer, 2015.
  * Rudin et al. [1992] Leonid Rudin, Stanley Osher, and Emad Fatemi.  Nonlinear total variation based noise removal algorithms.  _Physica D: Nonlinear Phenomena_ , 60:259–268, 11 1992.  doi: 10.1016/0167-2789(92)90242-F.
  * Saharia et al. [2022] Chitwan Saharia, William Chan, Huiwen Chang, Chris Lee, Jonathan Ho, Tim Salimans, David Fleet, and Mohammad Norouzi.  Palette: Image-to-image diffusion models.  In _ACM SIGGRAPH 2022 Conference Proceedings_ , pages 1–10, 2022.
  * Scaife [2020] AMM Scaife.  Big telescope, big data: towards exascale with the square kilometre array.  _Philosophical Transactions of the Royal Society A_ , 378(2166):20190060, 2020.
  * Schmidt et al. [2022] Kevin Schmidt, Felix Geyer, Stefan Fröse, P-S Blomenkamp, Marcus Brüggen, Francesco de Gasperin, Dominik Elsässer, and Wolfgang Rhode.  Deep learning-based imaging in radio interferometry.  _Astronomy & Astrophysics_, 664:A134, 2022.
  * Sohl-Dickstein et al. [2015] Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli.  Deep unsupervised learning using nonequilibrium thermodynamics.  In _International conference on machine learning_ , pages 2256–2265. PMLR, 2015.
  * Taran et al. [2023] O. Taran, O. Bait, M. Dessauges-Zavadsky, T. Holotyak, D. Schaerer, and S. Voloshynovskiy.  Challenging interferometric imaging: Machine learning-based source localization from uv-plane observations.  _Astronomy & Astrophysics_, may 2023.  doi: 10.1051/0004-6361/202245778.  URL https://doi.org/10.1051%2F0004-6361%2F202245778.
  * Tikhonov and Arsenin [1977] Andrey N. Tikhonov and Vasiliy Y. Arsenin.  _Solutions of ill-posed problems_.  V. H. Winston & Sons, Washington, D.C.: John Wiley & Sons, New York, 1977.  Translated from the Russian, Preface by translation editor Fritz John, Scripta Series in Mathematics.
  * Tingay et al. [2013] Steven John Tingay, Robert Goeke, Judd D Bowman, David Emrich, Stephen M Ord, Daniel A Mitchell, Miguel F Morales, Tom Booler, Brian Crosse, Randall B Wayth, et al.  The murchison widefield array: The square kilometre array precursor at low radio frequencies.  _Publications of the Astronomical Society of Australia_ , 30:e007, 2013.
  * Tolley et al. [2023] Emma Tolley, Simon Frasch, Etienne Orliac, Shreyam Krishna, Michele Bianco, Sepand Kashani, Paul Hurley, Matthieu Simeoni, and Jean-Paul Kneib.  Bipp: An efficient hpc implementation of the bluebild algorithm for radio astronomy, 2023.
  * Van der Walt et al. [2014] Stefan Van der Walt, Johannes L Schönberger, Juan Nunez-Iglesias, François Boulogne, Joshua D Warner, Neil Yager, Emmanuelle Gouillart, and Tony Yu.  scikit-image: image processing in python.  _PeerJ_ , 2:e453, 2014.
  * van Haarlem et al. [2013] Michael P van Haarlem, Michael W Wise, AW Gunst, George Heald, John P McKean, Jason WT Hessels, A Ger de Bruyn, Ronald Nijboer, John Swinbank, Richard Fallows, et al.  Lofar: The low-frequency array.  _Astronomy & astrophysics_, 556:A2, 2013.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.
  * Walter et al. [2016] Fabian Walter, Roberto Decarli, Manuel Aravena, Chris Carilli, Rychard Bouwens, Elisabete Da Cunha, Emanuele Daddi, Robert J Ivison, Dominik Riechers, Ian Smail, et al.  Alma spectroscopic survey in the hubble ultra deep field: survey description.  _The Astrophysical Journal_ , 833(1):67, 2016.
  * Wang et al. [2023] Ruoqi Wang, Zhuoyang Chen, Qiong Luo, and Feng Wang.  A conditional denoising diffusion probabilistic model for radio interferometric image reconstruction.  _arXiv preprint arXiv:2305.09121_ , 2023.
  * Wang et al. [2004] Zhou Wang, Alan C. Bovik, Hamid R. Sheikh, and Eero P. Simoncelli.  Image quality assessment: from error visibility to structural similarity.  _IEEE Transactions on Image Processing_ , 13(4):600–612, 2004.  URL http://dblp.uni-trier.de/db/journals/tip/tip13.html#WangBSS04.
  * Wiener [1949] Norbert Wiener.  _Extrapolation, interpolation, and smoothing of stationary time series, with engineering applications._ Technology Press of the Massachusetts Institute of Technology, Cambridge, 1949.

##  Appendix A Outliers

In Figures 14 and 15, we observe the outliers that persist irrespective of the
aggregation method used. We selected five sources exhibiting the greatest
relative error for further analysis. The problematic sources are those that
are located in close proximity to each other. The dataset includes two images
where two sources share identical coordinates, a scenario that the proposed
method cannot discern correctly. Instead of distinguishing the two, the model
predicts a single source with a flux equivalent to the sum of the fluxes from
the closely situated sources. The fifth error occurs when sources are located
near one another, but not with identical coordinates; while the flux of one is
accurately determined, the other is predicted to be substantially larger than
its actual value. This discrepancy can potentially be influenced by the
proximity of the neighboring source, as the predicted flux is of a comparable
magnitude. Detailed results are presented in Table 4.

Table 4: Most extreme outliers: sources with largest estimation errors in flux
by DDPM

Data Type | Ra | Dec | Flux ($\mu$Jy)
---|---|---|---
Real | 150.192944 | 3.995500 | 309
Real | 150.192944 | 3.995500 | 429
Predicted | 150.192945 | 3.995501 | 738
Real | 149.655167 | 3.643278 | 238
Real | 149.655167 | 3.643278 | 113
Predicted | 149.655165 | 3.643277 | 346
Real | 148.163222 | 1.370722 | 126
Real | 148.163222 | 1.370444 | 400
Predicted | 148.163220 | 1.370729 | 399
Predicted | 148.163223 | 1.370444 | 436

Note: These outliers are from the Figures 14 and 15. The proposed algorithm
fails when sources are overlapping each other. When there is exact overlap,
the prediction is a source with a flux that is approximately equal to the sum
of the fluxes of the overlapping true sources (first two sections).

##  Appendix B Diffusion without scaling

Figure 19: S/N versus reliability for $\gamma=1$.

This section presents plots for a DDPM model with $\gamma=1$ in preprocessing
during training in Equation 1. Figure 19 illustrates the absence of detected
low-S/N sources and a more uniformly spread reliability parameter compared to
$\gamma=2$.

Figure 20: Flux estimation for different aggregate-detect methods $\gamma=1$.
Figure 21: Flux estimation for detect-aggregate for $\gamma=1$.

Figures 20 and 21 show that the proposed model consistently underestimates
flux values for all ranges in a single run. Nevertheless, aggregation methods,
such as mean and median, significantly aid, unlike the detect-aggregate
method, which performs poorly, which is likely due to the consistent
underestimation of sources in a single run. Thus, directly averaging these
predictions does not improve the value. However, averaging over image pixels
(as aggregate-detect does) increases the source area, correcting overall flux
values.

Figure 22: Localization metrics for $\gamma=1$. Figure 23: Normalized metrics
for $\gamma=1$.

Figures 22 and 23 demonstrate the dependency of localization and
reconstruction metrics on the number of runs for $\gamma=1$. The
reconstruction metrics follow the same pattern as $\gamma=2$: a greater number
of images equates to better metrics. The localization metrics exhibit a
different behavior: a single run has very low completeness, less than 50%
compared to 96% for $\gamma=2$. For $\gamma=1$, the completes improves slower
than for $\gamma=2$, not surpassing 95% even after 20 runs, contrary to 98%
for $\gamma=2$.