Source code for ixpeobssim.evt.deconvolution

#!/usr/bin/env python
# Copyright (C) 2015--2020, the ixpeobssim team.
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

"""Module encapsulating image deconvolution facilities.

from __future__ import print_function, division

import numpy

from ixpeobssim.irf import DEFAULT_IRF_NAME, load_psf
from ixpeobssim.utils.logging_ import logger
from ixpeobssim.utils.units_ import degrees_to_arcmin

# pylint: disable=invalid-name

[docs] def circular_kernel(nside, normalize=False): """Build a simple circular kernel of a given side. Arguments --------- nside : int The side of the kernel array, assumed to be squared. normalize : bool If True, normalize to the kernel content (i.e., use for averaging rather than summing). """ shape = (nside, nside) kernel = numpy.zeros(shape, dtype=float) x, y = numpy.indices(shape) - nside / 2. + 0.5 kernel[x**2. + y**2. <= nside**2. / 4.] = 1. if normalize: kernel /= kernel.sum() return kernel
[docs] def calculate_psf_kernel(pixel_size, irf_name=DEFAULT_IRF_NAME, du_id=1, max_radius=0.03, sample_size=5000000): """Calculate the digital PSF kernel for image convolution. This is essentially creating a square grid whose pixel size is set to the value of the corresponding function argument and whose half size is (loosely) determined by the max_radius arguments. A number of points is thrown in the grid according to the PSF shape and the corresponding probability for each pixel is estimated by numerical integration. Arguments --------- pixel_size : float The pixels size of the image that the kernel is operating on (in decimal degrees). irf_name : str The IRF name. du_id : int The DU id. max_radius : float The maximum radius (in decimal degrees) for defining the size of the kernel grid. This should be large enough to fully contain the PSF. sample_size : int The size of the random sample for the Monte Carlo integration of the psf profile over the square grid. """'Building the PSF kernel for "%s" (DU %d)', irf_name, du_id) # Note the number of bins on both axes must be odd, so that the origin # of the PSF is at the center of a bin. nside = 2 * int(max_radius / pixel_size) + 1 delta = 0.5 * pixel_size * nside binning = numpy.linspace(-delta, delta, nside + 1) delta = degrees_to_arcmin(delta) args = nside, nside, delta'Sampling over a %d x %d spatial grid (+/-%.2f arcmin)', *args) # Load the PSF. psf = load_psf(irf_name, du_id) # Throw a whole nunch of random numbers with the PSF shape and bin them # over the kernel. delta = kernel, _ = numpy.histogramdd(delta, (binning, binning)) # Debug info to keep track of the kernel containment. kernel_sum = kernel.sum() containment = kernel.sum() / sample_size'Estimated PSF containment: %.9f', containment) # Normalize the kernel and we're done. return kernel / kernel_sum
def _kernel_radius(kernel): """Return the kernel radius in pixels. This is meant to return n, where the kernel shape is (2n + 1) x (2n + 1). """ nx, ny = kernel.shape assert nx == ny return nx // 2 def _kernel_mask(i, j, nk, image_shape): """Return the appropriate slice for convolving the generic pixel (i, j) within an image. For (i, j) = (50, 50) and a kernel radius of 2 pixels, e.g., this returns (slice(48, 53, None), slice(48, 53, None)). It is assumed that the kernel is square, with dimension (2nk + 1) x (2nk + 1). """ # Calculate the slice for the image---this is (i-nk--i+nk), (j-nk--j+nk) # modulo border effects. nx, ny = image_shape imin = max(i - nk, 0) imax = min(i + nk + 1, nx) jmin = max(j - nk, 0) jmax = min(j + nk + 1, ny) img_mask = (slice(imin, imax), slice(jmin, jmax)) # And now the slice for the kernel---this is (0--2nk+1), (0--2nk+1) # modulo border effects. nkmax = 2 * nk + 1 imin = max(0, nk - i) imax = min(nkmax, nx - i + nk) jmin = max(0, nk - j) jmax = min(nkmax, ny - j + nk) ker_mask = (slice(imin, imax), slice(jmin, jmax)) return img_mask, ker_mask
[docs] def convolve_image(image, kernel): """Convolve an image with the PSF kernel. While I am sure there exist a more clever (and faster) way to do this with scipy, I retained the explicit Python loop, here, as I don't think this is a bottleneck in any way, and having a full implementation is useful for debugging. Arguments --------- image : 2d array The input image to be convolved with the PSF. kernel : 2d array The PSF kernel for the convolution. """ output = numpy.full(image.shape, 0.) nk = _kernel_radius(kernel) for i, j in numpy.ndindex(image.shape): imask, kmask = _kernel_mask(i, j, nk, image.shape) output[imask] += kernel[kmask] * image[i, j] return output
[docs] def calculate_inverse_kernel(intensity, true_intensity, kernel): """Calculate the gargantuan (nx x ny x nk x nk) matrix that allows to deconvolve an (I, U or Q) image, given a proxy of the true intensity. The basic idea, here, is that if you have a sensible proxy of the true I map (e.g., a high-resolution image from Chandra) you can use it to gauge the fraction of the measured intensity ending up in each true pixel, and you can actually use that to deconvole a map (with the same binning) of either Q or U. Phrased in a slighlty different way, while the PSF kernel for the direct true -> measured convolution is identical for all pixels, the kernel for the inverse measured -> true deconvolution is in general different for all the pixels and depends on the true intensity. This is why we need to precompute a store a specific kernel for each pixels. This can be in turn use to deconvole any measured quantity with the same spatial binning (Q, U, W2). Warning ------- I am not positive this is right for the border pixels, where the kernel mask ends up out of the image---need to study that. Arguments --------- intensity : 2d array The measured intensity map. true_intensity : 2d array The proxy for the true intensity map. kernel : 2d array The PSF kernel for the convolution. """ # Before we start, we normalize the total intensity of the true image to that # of the reference image, since the two might, in principle, have totally # arbitrary normalizations. true_intensity *= intensity.sum() / true_intensity.sum() nk = _kernel_radius(kernel) K = numpy.full((*intensity.shape, *kernel.shape), 0.) for i, j in numpy.ndindex(intensity.shape): imask, kmask = _kernel_mask(i, j, nk, intensity.shape) subimg = intensity[imask] # Need to protect against zero-division errors. mask = (subimg != 0.) K[i, j][kmask][mask] = (kernel[kmask][mask] * true_intensity[i, j] / subimg[mask]) return K
[docs] def deconvolve_image(image, kernel, K): """Deconvolve a generic image (e.g., Q or U). Arguments --------- image : 2d array The input image to be deconvoled. kernel : 2d array The PSF kernel. K : 4d array The inverse pixel-by-pixel kernel for the deconvolution. """ output = numpy.full(image.shape, 0.) nk = _kernel_radius(kernel) for i, j in numpy.ndindex(image.shape): imask, kmask = _kernel_mask(i, j, nk, image.shape) output[i, j] = (image[imask] * K[i, j][kmask]).sum() return output