**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.