Authors Aneta Siemiginowska Axel Donath David van Dyk Douglas Burke Karthik Reddy Solipuram Vinay Kashyap
License CC-BY-4.0
Creative Commons Attribution 4.0 International (CC BY 4.0)
https://creativecommons.org/licenses/by/4.0/
Access to this work was provided by the University of Maryland, Baltimore County (UMBC)
ScholarWorks@UMBC digital repository on the Maryland Shared Open Access (MD-SOAR)
platform.
Please provide feedback
Please support the ScholarWorks@UMBC repository by emailing scholarworks-
group@umbc.edu and telling us what having access to this work means to you and why it’s
important to you. Thank you.
98 PROC. OF THE 21st PYTHON IN SCIENCE CONF. (SCIPY 2022)
Pylira: deconvolution of images in the presence of
Poisson noise
Axel Donath‡∗ , Aneta Siemiginowska‡ , Vinay Kashyap‡ , Douglas Burke‡ , Karthik Reddy Solipuram§ , David van Dyk¶
F
Abstract—All physical and astronomical imaging observations are degraded by of the signal intensity to the signal variance. Any statistically
the finite angular resolution of the camera and telescope systems. The recovery correct post-processing or reconstruction method thus requires a
of the true image is limited by both how well the instrument characteristics careful treatment of the Poisson nature of the measured image.
are known and by the magnitude of measurement noise. In the case of a To maximise the scientific use of the data, it is often desired to
high signal to noise ratio data, the image can be sharpened or “deconvolved”
correct the degradation introduced by the imaging process. Besides
robustly by using established standard methods such as the Richardson-Lucy
method. However, the situation changes for sparse data and the low signal to
correction for non-uniform exposure and background noise this
noise regime, such as those frequently encountered in X-ray and gamma-ray also includes the correction for the "blurring" introduced by the
astronomy, where deconvolution leads inevitably to an amplification of noise point spread function (PSF) of the instrument. Where the latter
and poorly reconstructed images. However, the results in this regime can process is often called "deconvolution". Depending on whether
be improved by making use of physically meaningful prior assumptions and the PSF of the instrument is known or not, one distinguishes
statistically principled modeling techniques. One proposed method is the LIRA between the "blind deconvolution" and "non blind deconvolution"
algorithm, which requires smoothness of the reconstructed image at multiple process. For astronomical observations, the PSF can often either
scales. In this contribution, we introduce a new python package called Pylira,
be simulated, given a model of the telescope and detector, or
which exposes the original C implementation of the LIRA algorithm to Python
inferred directly from the data by observing far distant objects,
users. We briefly describe the package structure, development setup and show
a Chandra as well as Fermi-LAT analysis example.
which appear as a point source to the instrument.
While in other branches of astronomy deconvolution methods
Index Terms—deconvolution, point spread function, poisson, low counts, X-ray, are already part of the standard analysis, such as the CLEAN
gamma-ray algorithm for radio data, developed by [Hog74], this is not the
case for X-ray and gamma-ray astronomy. As any deconvolution
method aims to enhance small-scale structures in an image, it
Introduction becomes increasingly hard to solve for the regime of low signal-
Any physical and astronomical imaging process is affected by to-noise ratio, where small-scale structures are more affected by
the limited angular resolution of the instrument or telescope. In noise.
addition, the quality of the resulting image is also degraded by
background or instrumental measurement noise and non-uniform The Deconvolution Problem
exposure. For short wavelengths and associated low intensities of Basic Statistical Model
the signal, the imaging process consists of recording individual Assuming the data in each pixel di in the recorded counts image
photons (often called "events") originating from a source of follows a Poisson distribution, the total likelihood of obtaining the
interest. This imaging process is typical for X-ray and gamma- measured image from a model image of the expected counts λi
ray telescopes, but images taken by magnetic resonance imaging with N pixels is given by:
or fluorescence microscopy show Poisson noise too. For each
individual photon, the incident direction, energy and arrival time
N exp −di λidi
L (d|λ ) = ∏ (1)
is measured. Based on this information, the event can be binned i di !
into two dimensional data structures to form an actual image.
By taking the logarithm, dropping the constant terms and inverting
As a consequence of the low intensities associated to the the sign one can transform the product into a sum over pixels,
recording of individual events, the measured signal follows Pois- which is also often called the Cash [Cas79] fit statistics:
son statistics. This imposes a non-linear relationship between the
N
measured signal and true underlying intensity as well as a coupling C (λ |d) = ∑(λi − di log λi ) (2)
i
* Corresponding author: axel.donath@cfa.harvard.edu
‡ Center for Astrophysics | Harvard & Smithsonian Where the expected counts λi are given by the convolution of the
§ University of Maryland Baltimore County true underlying flux distribution xi with the PSF pk :
¶ Imperial College London
λi = ∑ xi pi−k (3)
Copyright © 2022 Axel Donath et al. This is an open-access article distributed k
under the terms of the Creative Commons Attribution License, which permits
unrestricted use, distribution, and reproduction in any medium, provided the This operation is often called "forward modelling" or "forward
original author and source are credited. folding" with the instrument response.
PYLIRA: DECONVOLUTION OF IMAGES IN THE PRESENCE OF POISSON NOISE 99
Richardson Lucy (RL)
To obtain the most likely value of xn given the data, one searches
a maximum of the total likelihood function, or equivalently a of
minimum C . This high dimensional optimization problem can
e.g., be solved by a classic gradient descent approach. Assuming
the pixels values xi of the true image as independent parameters,
one can take the derivative of Eq. 2 with respect to the individual
xi . This way one obtains a rule for how to update the current set
of pixels xn in each iteration of the optimization:
∂ C (d|x)
xn+1 = xn − α · (4)
∂ xi
Where α is a factor to define the step size. This method is in
general equivalent to the gradient descent and backpropagation
methods used in modern machine learning techniques. This ba-
sic principle of solving the deconvolution problem for images
with Poisson noise was proposed by [Ric72] and [Luc74]. Their
method, named after the original authors, is often known as the
Fig. 1: The images show the result of the RL algorithm applied
Richardson & Lucy (RL) method. It was shown by [Ric72] that to a simulated example dataset with varying numbers of iterations.
this converges to a maximum likelihood solution of Eq. 2. A The image in the upper left shows the simulated counts. Those have
Python implementation of the standard RL method is available been derived from the ground truth (upper mid) by convolving with a
e.g. in the Scikit-Image package [vdWSN+ 14]. Gaussian PSF of width σ = 3 pix and applying Poisson noise to it.
Instead of the iterative, gradient descent based optimization it The illustration uses the implementation of the RL algorithm from the
is also possible to sample from the posterior distribution using a Scikit-Image package [vdWSN+ 14].
simple Metropolis-Hastings [Has70] approach and uniform prior.
This is demonstrated in one of the Pylira online tutorials (Intro-
the smoothness of the reconstructed image on multiple spatial
duction to Deconvolution using MCMC Methods).
scales. Starting from the full resolution, the image pixels xi are
collected into 2 by 2 groups Qk . The four pixel values associated
RL Reconstruction Quality
with each group are divided by their sum to obtain a grid of “split
While technically the RL method converges to a maximum like- proportions” with respect to the image down-sized by a factor of
lihood solution, it mostly still results in poorly restored images, two along both axes. This process is repeated using the down sized
especially if extended emission regions are present in the image. image with pixel values equal to the sums over the 2 by 2 groups
The problem is illustrated in Fig. 1 using a simulated example from the full-resolution image, and the process continues until the
image. While for a low number of iterations, the RL method still resolution of the image is only a single pixel, containing the total
results in a smooth intensity distribution, the structure of the image sum of the full-resolution image. This multi-scale representation
decomposes more and more into a set of point-like sources with is illustrated in Fig. 2.
growing number of iterations. For each of the 2x2 groups of the re-normalized images a
Because of the PSF convolution, an extended emission region Dirichlet distribution is introduced as a prior:
can decompose into multiple nearby point sources and still lead
to good model prediction, when compared with the data. Those φk ∝ Dirichlet(αk , αk , αk , αk ) (6)
almost equally good solutions correspond to many narrow local and multiplied across all 2x2 groups and resolution levels k. For
minima or "spikes" in the global likelihood surface. Depending on each resolution level a smoothing parameter αk is introduced.
the start estimate for the reconstructed image x the RL method These hyper-parameters can be interpreted as having an infor-
will follow the steepest gradient and converge towards the nearest mation content equivalent of adding αk "hallucinated" counts in
narrow local minimum. This problem has been described by each grouping. This effectively results in a smoothing of the
multiple authors, such as [PR94] and [FBPW95]. image at the given resolution level. The distribution of α values
at each resolution level is the further described by a hyper-prior
Multi-Scale Prior & LIRA
distribution:
One solution to this problem was described in [ECKvD04] and p(αk ) = exp (−δ α 3 /3) (7)
[CSv+ 11]. First, the simple forward folded model described in
Eq. 3 can be extended by taking into account the non-uniform Resulting in a fully hierarchical Bayesian model. A more com-
exposure ei and an additional known background component bi : plete and detailed description of the prior definition is given in
[ECKvD04].
λi = ∑ (ei · (xi + bi )) pi−k (5) The problem is then solved by using a Gibbs MCMC sampling
k approach. After a "burn-in" phase the sampling process typically
The background bi can be more generally understood as a "base- reaches convergence and starts sampling from the posterior distri-
line" image and thus include known structures, which are not of bution. The reconstructed image is then computed as the mean of
interest for the deconvolution process. E.g., a bright point source the posterior samples. As for each pixel a full distribution of its
to model the core of an AGN while studying its jets. values is available, the information can also be used to compute
Second, the authors proposed to extend the Poisson log- the associated error of the reconstructed value. This is another
likelihood function (Equation 2) by a log-prior term that controls main advantage over RL or Maxium A-Postori (MAP) algorithms.
100 PROC. OF THE 21st PYTHON IN SCIENCE CONF. (SCIPY 2022)
1 $ sudo apt-get install r-base-dev r-base r-mathlib
2 $ pip install pylira
For more detailed instructions see Pylira installation instructions.
API & Subpackages
Pylira is structured in multiple sub-packages. The pylira.src
module contains the original C implementation and the Pybind11
wrapper code. The pylira.core sub-package contains the
main Python API, pylira.utils includes utility functions
for plotting and serialisation. And pylira.data implements
multiple pre-defined datasets for testing and tutorials.
Analysis Examples
Simple Point Source
Pylira was designed to offer a simple Python class based user
interface, which allows for a short learning curve of using the
package for users who are familiar with Python in general and
more specifically with Numpy. A typical complete usage example
of the Pylira package is shown in the following:
Fig. 2: The image illustrates the multi-scale decomposition used in
the LIRA prior for a 4x4 pixels example image. Each quadrant of 2x2 1 import numpy as np
sub-images is labelled with QN . The sub-pixels in each quadrant are 2 from pylira import LIRADeconvolver
labelled Λi j . . 3 from pylira.data import point_source_gauss_psf
4
5 # create example dataset
6 data = point_source_gauss_psf()
The Pylira Package 7
8 # define initial flux image
Dependencies & Development 9 data["flux_init"] = data["flux"]
The Pylira package is a thin Python wrapper around the original 10
LIRA implementation provided by the authors of [CSv+ 11]. The 11 deconvolve = LIRADeconvolver(
12 n_iter_max=3_000,
original algorithm was implemented in C and made available as a 13 n_burn_in=500,
package for the R Language [R C20]. Thus the implementation de- 14 alpha_init=np.ones(5)
pends on the RMath library, which is still a required dependency of 15 )
16
Pylira. The Python wrapper was built using the Pybind11 [JRM17] 17 result = deconvolve.run(data=data)
package, which allows to reduce the code overhead introduced by 18
the wrapper to a minimum. For the data handling, Pylira relies on 19 # plot pixel traces, result shown in Figure 3
Numpy [HMvdW+ 20] arrays for the serialisation to the FITS data 20 result.plot_pixel_traces_region(
21 center_pix=(16, 16), radius_pix=3
format on Astropy [Col18]. The (interactive) plotting functionality 22 )
is achieved via Matplotlib [Hun07] and Ipywidgets [wc15], which 23
are both optional dependencies. Pylira is openly developed on 24 # plot pixel traces, result shown in Figure 4
25 result.plot_parameter_traces()
Github at https://github.com/astrostat/pylira. It relies on GitHub 26
Actions as a continuous integration service and uses the Read 27 # finally serialise the result
the Docs service to build and deploy the documentation. The on- 28 result.write("result.fits")
line documentation can be found on https://pylira.readthedocs.io. The main interface is exposed via the LIRADeconvolver
Pylira implements a set of unit tests to assure compatibility class, which takes the configuration of the algorithm on initial-
and reproducibility of the results with different versions of the isation. Typical configuration parameters include the total num-
dependencies and across different platforms. As Pylira relies on ber of iterations n_iter_max and the number of "burn-in"
random sampling for the MCMC process an exact reproducibility iterations, to be excluded from the posterior mean computation.
of results is hard to achieve on different platforms; however the The data, represented by a simple Python dict data structure,
agreement of results is at least guaranteed in the statistical limit of contains a "counts", "psf" and optionally "exposure"
drawing many samples. and "background" array. The dataset is then passed to the
LIRADeconvolver.run() method to execute the deconvolu-
Installation
tion. The result is a LIRADeconvolverResult object, which
Pylira is available via the Python package index (pypi.org), features the possibility to write the result as a FITS file, as well
currently at version 0.1. As Pylira still depends on the RMath as to inspect the result with diagnostic plots. The result of the
library, it is required to install this first. So the recommended way computation is shown in the left panel of Fig. 3.
to install Pylira is on MacOS is:
1 $ brew install r Diagnostic Plots
2 $ pip install pylira
To validate the quality of the results Pylira provides many built-
On Linux the RMath dependency can be installed using standard in diagnostic plots. One of these diagnostic plot is shown in the
package managers. For example on Ubuntu, one would do right panel of Fig. 3. The plot shows the image sampling trace
PYLIRA: DECONVOLUTION OF IMAGES IN THE PRESENCE OF POISSON NOISE 101
Pixel trace for (16, 16)
30 800 1000
700
25 800
600
20
500 600
Posterior Mean
Burn in
Valid
15 400 Mean
400 1 Std. Deviation
300
10
200 200
5
100
0
0
0 5 10 15 20 25 30 0 500 1000 1500 2000 2500 3000
Number of Iterations
Fig. 3: The curves show the traces of value the pixel of interest for a simulated point source and its neighboring pixels (see code example).
The image on the left shows the posterior mean. The white circle in the image shows the circular region defining the neighboring pixels. The
blue line on the right plot shows the trace of the pixel of interest. The solid horizontal orange line shows the mean value (excluding burn-in)
of the pixel across all iterations and the shaded orange area the 1 σ error region. The burn in phase is shown in transparent blue and ignored
while computing the mean. The shaded gray lines show the traces of the neighboring pixels.
for a single pixel of interest and its surrounding circular region of Chandra is a space-based X-ray observatory, which has been
interest. This visualisation allows the user to assess the stability in operation since 1999. It consists of nested cylindrical paraboloid
of a small region in the image e.g. an astronomical point source and hyperboloid surfaces, which form an imaging optical system
during the MCMC sampling process. Due to the correlation with for X-rays. In the focal plane, it has multiple instruments for dif-
neighbouring pixels, the actual value of a pixel might vary in the ferent scientific purposes. This includes a high-resolution camera
sampling process, which appears as "dips" in the trace of the pixel (HRC) and an Advanced CCD Imaging Spectrometer (ACIS). The
of interest and anti-correlated "peaks" in the one or mutiple of typical angular resolution is 0.5 arcsecond and the covered energy
the surrounding pixels. In the example a stable state of the pixels ranges from 0.1 - 10 keV.
of interest is reached after approximately 1000 iterations. This Figure 5 shows the result of the Pylira algorithm applied to
suggests that the number of burn-in iterations, which was defined Chandra data of the Galactic Center region between 0.5 and 7 keV.
beforehand, should be increased. The PSF was obtained from simulations using the simulate_psf
Pylira relies on an MCMC sampling approach to sample tool from the official Chandra science tools ciao 4.14 [FMA+ 06].
a series of reconstructed images from the posterior likelihood The algorithm achieves both an improved spatial resolution as well
defined by Eq. 2. Along with the sampling, it marginalises over as a reduced noise level and higher contrast of the image in the
the smoothing hyper-parameters and optimizes them in the same right panel compared to the unprocessed counts data shown in the
process. To diagnose the validity of the results it is important to left panel.
visualise the sampling traces of both the sampled images as well As a second example, we use data from the Fermi Large Area
as hyper-parameters. Telescope (LAT). The Fermi-LAT is a satellite-based imaging
Figure 4 shows another typical diagnostic plot created by the gamma-ray detector, which covers an energy range of 20 MeV
code example above. In a multi-panel figure, the user can inspect to >300 GeV. The angular resolution varies strongly with energy
the traces of the total log-posterior as well as the traces of the and ranges from 0.1 to >10 degree1 .
smoothing parameters. Each panel corresponds to the smoothing Figure 6 shows the result of the Pylira algorithm applied to
hyper parameter introduced for each level of the multi-scale Fermi-LAT data above 1 GeV to the region around the Galactic
representation of the reconstructed image. The figure also shows Center. The PSF was obtained from simulations using the gtpsf
the mean value along with the 1 σ error region. In this case, tool from the official Fermitools v2.0.19 [Fer19]. First, one can
the algorithm shows stable convergence after a burn-in phase of see that the algorithm achieves again a considerable improvement
approximately 200 iterations for the log-posterior as well as all of in the spatial resolution compared to the raw counts. It clearly
the multi-scale smoothing parameters. resolves multiple point sources left to the bright Galactic Center
source.
Astronomical Analysis Examples Summary & Outlook
Both in the X-ray as well as in the gamma-ray regime, the Galactic The Pylira package provides Python wrappers for the LIRA al-
Center is a complex emission region. It shows point sources, gorithm. It allows the deconvolution of low-counts data following
extended sources, as well as underlying diffuse emission and thus 1. https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.
represents a challenge for any astronomical data analysis. htm
102 PROC. OF THE 21st PYTHON IN SCIENCE CONF. (SCIPY 2022)
Logpost Smoothingparam0 Smoothingparam1
Burn in 0.35 0.35
1500 Valid 0.30
Mean 0.30
1 Std. Deviation 0.25 0.25
1000
0.20 0.20
500 0.15 0.15
0 0.10 0.10
0.05 0.05
500
0.00 0.00
0 200 400 600 800 1000 0 200 400 600 800 1000 0 200 400 600 800 1000
Number of Iterations Number of Iterations Number of Iterations
Smoothingparam2 Smoothingparam3 Smoothingparam4
0.200
0.175
0.20 0.175
0.150 0.150
0.15 0.125 0.125
0.100 0.100
0.10 0.075
0.075
0.05 0.050 0.050
0.025 0.025
0.00 0.000 0.000
0 200 400 600 800 1000 0 200 400 600 800 1000 0 200 400 600 800 1000
Number of Iterations Number of Iterations Number of Iterations
Fig. 4: The curves show the traces of the log posterior value as well as traces of the values of the prior parameter values. The SmoothingparamN
parameters correspond to the smoothing parameters αN per multi-scale level. The solid horizontal orange lines show the mean value, the shaded
orange area the 1 σ error region. The burn in phase is shown transparent and ignored while estimating the mean.
Counts Deconvolved
500
PSF
257
132
-29°00'25"
68
Declination
35
Counts
18
30"
9
5
2
35"
17h45m40.6s40.4s 40.2s 40.0s 39.8s 39.6s 17h45m40.6s40.4s 40.2s 40.0s 39.8s 39.6s
Right Ascension Right Ascension
Fig. 5: Pylira applied to Chandra ACIS data of the Galactic Center region, using the observation IDs 4684 and 4684. The image on the left
shows the raw observed counts between 0.5 and 7 keV. The image on the right shows the deconvolved version. The LIRA hyperprior values
were chosen as ms_al_kap1=1, ms_al_kap2=0.02, ms_al_kap3=1. No baseline background model was included.
PYLIRA: DECONVOLUTION OF IMAGES IN THE PRESENCE OF POISSON NOISE 103
Counts Deconvolved
200
0°40' PSF
120
72
20'
43
Galactic Latitude
00' 26
Counts
16
-0°20' 9
5
40' 2
0°40' 20' 00' 359°40' 20' 0°40' 20' 00' 359°40' 20'
Galactic Longitude Galactic Longitude
Fig. 6: Pylira applied to Fermi-LAT data from the Galactic Center region. The image on the left shows the raw measured counts between
5 and 1000 GeV. The image on the right shows the deconvolved version. The LIRA hyperprior values were chosen as ms_al_kap1=1,
ms_al_kap2=0.02, ms_al_kap3=1. No baseline background model was included.
Poisson statistics using a Bayesian sampling approach and a multi- [CSv+ 11] A. Connors, N. M. Stein, D. van Dyk, V. Kashyap, and
scale smoothing prior assumption. The results can be easily written A. Siemiginowska. LIRA — The Low-Counts Image Restora-
tion and Analysis Package: A Teaching Version via R. In I. N.
to FITS files and inspected by plotting the trace of the sampling Evans, A. Accomazzi, D. J. Mink, and A. H. Rots, editors,
process. This allows users to check for general convergence as Astronomical Data Analysis Software and Systems XX, volume
well as pixel to pixel correlations for selected regions of interest. 442 of Astronomical Society of the Pacific Conference Series,
The package is openly developed on GitHub and includes tests page 463, July 2011.
[ECKvD04] David N. Esch, Alanna Connors, Margarita Karovska, and
and documentation, such that it can be maintained and improved David A. van Dyk. An image restoration technique with
in the future, while ensuring consistency of the results. It comes error estimates. The Astrophysical Journal, 610(2):1213–
with multiple built-in test datasets and explanatory tutorials in 1227, aug 2004. URL: https://doi.org/10.1086/421761, doi:
10.1086/421761.
the form of Jupyter notebooks. Future plans include the support [FBPW95] D. A. Fish, A. M. Brinicombe, E. R. Pike, and J. G.
for parallelisation or distributed computing, more flexible prior Walker. Blind deconvolution by means of the richardson–
definitions and the possibility to account for systematic errors on lucy algorithm. J. Opt. Soc. Am. A, 12(1):58–65, Jan 1995.
the PSF during the sampling process. URL: http://opg.optica.org/josaa/abstract.cfm?URI=josaa-12-
1-58, doi:10.1364/JOSAA.12.000058.
[Fer19] Fermi Science Support Development Team. Fermitools: Fermi
Acknowledgements Science Tools. Astrophysics Source Code Library, record
ascl:1905.011, May 2019. arXiv:1905.011.
This work was conducted under the auspices of the CHASC [FMA+ 06] Antonella Fruscione, Jonathan C. McDowell, Glenn E. Allen,
International Astrostatistics Center. CHASC is supported by NSF Nancy S. Brickhouse, Douglas J. Burke, John E. Davis, Nick
Durham, Martin Elvis, Elizabeth C. Galle, Daniel E. Har-
grants DMS-21-13615, DMS-21-13397, and DMS-21-13605; by ris, David P. Huenemoerder, John C. Houck, Bish Ishibashi,
the UK Engineering and Physical Sciences Research Council Margarita Karovska, Fabrizio Nicastro, Michael S. Noble,
[EP/W015080/1]; and by NASA 18-APRA18-0019. We thank Michael A. Nowak, Frank A. Primini, Aneta Siemiginowska,
CHASC members for many helpful discussions, especially Xiao- Randall K. Smith, and Michael Wise. CIAO: Chandra’s data
analysis system. In David R. Silva and Rodger E. Doxsey,
Li Meng and Katy McKeough. DvD was also supported in part editors, Society of Photo-Optical Instrumentation Engineers
by a Marie-Skodowska-Curie RISE Grant (H2020-MSCA-RISE- (SPIE) Conference Series, volume 6270 of Society of Photo-
2019-873089) provided by the European Commission. Aneta Optical Instrumentation Engineers (SPIE) Conference Series,
page 62701V, June 2006. doi:10.1117/12.671760.
Siemiginowska, Vinay Kashyap, and Doug Burke further acknowl- [Has70] W. K. Hastings. Monte Carlo Sampling Methods using Markov
edge support from NASA contract to the Chandra X-ray Center Chains and their Applications. Biometrika, 57(1):97–109,
NAS8-03060. April 1970. doi:10.1093/biomet/57.1.97.
[HMvdW+ 20] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der
Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric
R EFERENCES Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith,
Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van
[Cas79] W. Cash. Parameter estimation in astronomy through ap- Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del
plication of the likelihood ratio. The Astrophysical Journal, Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant,
228:939–947, March 1979. doi:10.1086/156922. Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer
[Col18] Astropy Collaboration. The Astropy Project: Building an Abbasi, Christoph Gohlke, and Travis E. Oliphant. Array pro-
Open-science Project and Status of the v2.0 Core Package. The gramming with NumPy. Nature, 585(7825):357–362, Septem-
Astrophysical Journal, 156(3):123, September 2018. arXiv: ber 2020. URL: https://doi.org/10.1038/s41586-020-2649-2,
1801.02634, doi:10.3847/1538-3881/aabc4f. doi:10.1038/s41586-020-2649-2.
104 PROC. OF THE 21st PYTHON IN SCIENCE CONF. (SCIPY 2022)
[Hog74] J. A. Hogbom. Aperture Synthesis with a Non-Regular
Distribution of Interferometer Baselines. Astronomy and As-
trophysics Supplement, 15:417, June 1974.
[Hun07] J. D. Hunter. Matplotlib: A 2d graphics environment. Com-
puting in Science & Engineering, 9(3):90–95, 2007. doi:
10.1109/MCSE.2007.55.
[JRM17] Wenzel Jakob, Jason Rhinelander, and Dean Moldovan. py-
bind11 – seamless operability between c++11 and python,
2017. https://github.com/pybind/pybind11.
[Luc74] L. B. Lucy. An iterative technique for the rectification of
observed distributions. Astronomical Journal, 79:745, June
1974. doi:10.1086/111605.
[PR94] K. M. Perry and S. J. Reeves. Generalized Cross-Validation
as a Stopping Rule for the Richardson-Lucy Algorithm. In
Robert J. Hanisch and Richard L. White, editors, The Restora-
tion of HST Images and Spectra - II, page 97, January 1994.
doi:10.1002/ima.1850060412.
[R C20] R Core Team. R: A Language and Environment for Statistical
Computing. R Foundation for Statistical Computing, Vienna,
Austria, 2020. URL: https://www.R-project.org/.
[Ric72] William Hadley Richardson. Bayesian-Based Iterative Method
of Image Restoration. Journal of the Optical Society of
America (1917-1983), 62(1):55, January 1972. doi:10.
1364/josa.62.000055.
[vdWSN+ 14] Stéfan van der Walt, Johannes L. Schönberger, Juan Nunez-
Iglesias, François Boulogne, Joshua D. Warner, Neil Yager,
Emmanuelle Gouillart, Tony Yu, and the scikit-image con-
tributors. scikit-image: image processing in Python. PeerJ,
2:e453, 6 2014. URL: https://doi.org/10.7717/peerj.453, doi:
10.7717/peerj.453.
[wc15] Jupyter widgets community. ipywidgets, a github repository.
Retrieved from https://github.com/jupyter-widgets/ipywidgets,
2015.