DOKK Library

Pylira: deconvolution of images in the presence of Poisson noise

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)

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)

Please provide feedback
Please support the ScholarWorks@UMBC repository by emailing scholarworks- 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¶


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
measured signal and true underlying intensity as well as a coupling                                       C (λ |d) = ∑(λi − di log λi )                   (2)
* Corresponding author:
‡ 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
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
                     λ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 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 import point_source_gauss_psf
                                                                           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   )
Pylira. The Python wrapper was built using the Pybind11 [JRM17]           17   result =
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 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         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
                                                                 method to execute the deconvolu-
                                                                          tion. The result is a LIRADeconvolverResult object, which
Pylira is available via the Python package index (,              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
    25                                                                            800
                                                              500                 600

                                                                 Posterior Mean
                                                                                                                                           Burn in
    15                                                        400                                                                          Mean
                                                                                  400                                                      1 Std. Deviation
                                                              200                 200
         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

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.
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
                                                          0.20                                        0.20
           500                                            0.15                                        0.15
                0                                         0.10                                        0.10
                                                          0.05                                        0.05
                                                          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.20                                            0.175
                                                         0.150                                       0.150
         0.15                                            0.125                                       0.125
                                                         0.100                                       0.100
         0.10                                                                                        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



                        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
                          0°40'            PSF
      Galactic Latitude

                            00'                                                                                                                           26


                          -0°20'                                                                                                                          9
                            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:, doi:
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:
                                                                                                         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:,
                             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:
[JRM17]     Wenzel Jakob, Jason Rhinelander, and Dean Moldovan. py-
            bind11 – seamless operability between c++11 and python,
[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.
[R C20]     R Core Team. R: A Language and Environment for Statistical
            Computing. R Foundation for Statistical Computing, Vienna,
            Austria, 2020. URL:
[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.
[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:, doi:
[wc15]      Jupyter widgets community. ipywidgets, a github repository.
            Retrieved from,