# Licensed under a 3-clause BSD style license - see LICENSE.rst
import copy
import logging
import numpy as np
from scipy.signal import fftconvolve
from scipy.ndimage import label
from astropy.io import fits
from astropy.table import Table
from astropy.convolution import Gaussian2DKernel, MexicanHat2DKernel
from ..maps import WcsNDMap, MapAxis, WcsGeom
__all__ = ["CWT", "CWTData", "CWTKernels"]
log = logging.getLogger(__name__)
def difference_of_gauss_kernel(radius, scale_step, n_sigmas=8):
"""Difference of 2 Gaussians (i.e. Mexican hat) kernel array.
TODO: replace by http://astropy.readthedocs.io/en/latest/api/astropy.convolution.MexicanHat2DKernel.html
once there are tests in place that establish the algorithm
"""
sizex = int(n_sigmas * scale_step * radius)
sizey = int(n_sigmas * scale_step * radius)
radius = float(radius)
xc = 0.5 * sizex
yc = 0.5 * sizey
y, x = np.mgrid[0 : sizey - 1, 0 : sizex - 1]
x = x - xc
y = y - yc
x1 = x / radius
y1 = y / radius
g1 = np.exp(-0.5 * (x1 ** 2 + y1 ** 2))
g1 = g1 / (2 * np.pi * radius ** 2)
x1 = x1 / scale_step
y1 = y1 / scale_step
g2 = np.exp(-0.5 * (x1 ** 2 + y1 ** 2))
g2 = g2 / (2 * np.pi * radius ** 2 * scale_step ** 2)
return g1 - g2
[docs]class CWT:
"""Continuous wavelet transform.
TODO: describe algorithm (modify the words below)
Depending on their spectral index, sources won't have the same characteristic scale.
Therefore to detect sources, we need to compute the wavelet transform at several
scales in order to search for various PSF sizes. Then for each scale, the wavelet
transform values under the given significance threshold are rejected. This gives us
a multiscale support. Then, using the reconstruction by continuous wavelet packets,
we obtain a filtered image yielding the detected sources.
To compute the threshold image for a given scale a, the standard EGRET diffuse
background model to which was added the flux of the extragalactic background,
and the exposure for the considered energy range were used.
Parameters
----------
kernels : `~gammapy.detect.CWTKernels`
Kernels for the algorithm.
max_iter : int, optional (default 10)
The maximum number of iterations of the CWT algorithm.
tol : float, optional (default 1e-5)
Tolerance for stopping criterion.
significance_threshold : float, optional (default 3.0)
Measure of statistical significance.
significance_island_threshold : float, optional (default None)
Measure is used for cleaning of isolated pixel islands
that are slightly above ``significance_threshold``.
remove_isolated : boolean, optional (default True)
If ``True``, isolated pixels will be removed.
keep_history : boolean, optional (default False)
Save cwt data from all the iterations.
References
----------
R. Terrier et al (2001) "Wavelet analysis of EGRET data"
See https://ui.adsabs.harvard.edu/abs/2001ICRC....7.2923T
"""
def __init__(
self,
kernels,
max_iter=10,
tol=1e-5,
significance_threshold=3.0,
significance_island_threshold=None,
remove_isolated=True,
keep_history=False,
):
self.kernels = kernels
self.max_iter = max_iter
self.tol = tol
self.significance_threshold = significance_threshold
self.significance_island_threshold = significance_island_threshold
self.remove_isolated = remove_isolated
self.history = [] if keep_history else None
# previous_variance is initialized on the first iteration
self.previous_variance = None
def _execute_iteration(self, data):
"""Do one iteration of the algorithm.
Parameters
----------
data : `~gammapy.detect.CWTData`
Images.
"""
self._transform(data=data)
self._compute_support(data=data)
self._inverse_transform(data=data)
def _transform(self, data):
"""Do the transform itself.
The transform is made by using `scipy.signal.fftconvolve`.
TODO: document.
Parameters
----------
data : `~gammapy.detect.CWTData`
Images for transform.
"""
total_background = data._model + data._background + data._approx
excess = data._counts - total_background
log.debug("Excess sum: {:.4f}".format(excess.sum()))
log.debug("Excess max: {:.4f}".format(excess.max()))
log.debug("Computing transform and error")
for idx_scale, kern in self.kernels.kern_base.items():
data._transform_3d[idx_scale] = fftconvolve(excess, kern, mode="same")
data._error[idx_scale] = np.sqrt(
fftconvolve(total_background, kern ** 2, mode="same")
)
log.debug("Error sum: {:.4f}".format(data._error.sum()))
log.debug("Error max: {:.4f}".format(data._error.max()))
log.debug("Computing approx and approx_bkg")
data._approx = fftconvolve(
data._counts - data._model - data._background,
self.kernels.kern_approx,
mode="same",
)
data._approx_bkg = fftconvolve(
data._background, self.kernels.kern_approx, mode="same"
)
log.debug("Approximate sum: {:.4f}".format(data._approx.sum()))
log.debug("Approximate background sum: {:.4f}".format(data._approx_bkg.sum()))
def _compute_support(self, data):
"""Compute the multiresolution support with hard sigma clipping.
Imposing a minimum significance on a connected region of significant pixels
(i.e. source detection).
TODO: document?
What's happening here:
- calc significance for all scales
- for each scale compute support:
- create mask, whether significance value more than significance_threshold
- find all separated structures on mask, label them (with help
of 'scipy.ndimage.label')
- for each structure do:
- find pixels coordinates of the structure
- if just one pixel, we can remove it from mask if we want
- if the max value of significance in that structure less
than significance_island_threshold, remove this structure from mask
- update support by mask
Parameters
----------
data : `~gammapy.detect.CWTData`
Images after transform.
"""
log.debug("Computing significance")
significance = data._transform_3d / data._error
log.debug("For each scale start to compute support")
for idx_scale in range(self.kernels.n_scale):
log.debug(
"Start to compute support for scale "
"{:.2f}".format(self.kernels.scales[idx_scale])
)
log.debug(
"Create mask based on significance "
"threshold {:.2f}".format(self.significance_threshold)
)
mask = significance[idx_scale] > self.significance_threshold
# Produce a list of connex structures in the support
labeled_mask, n_structures = label(mask)
for struct_label in range(n_structures):
coords = np.where(labeled_mask == struct_label + 1)
# Remove isolated pixels from support
if self.remove_isolated and coords[0].size == 1:
log.debug("Remove isolated pixels from support")
mask[coords] = False
if self.significance_island_threshold is not None:
# If maximal significance of the structure does not reach significance
# island threshold, remove significant pixels island from support
struct_signif = significance[idx_scale][coords]
if struct_signif.max() < self.significance_island_threshold:
log.debug(
"Remove significant pixels island {} from support".format(
struct_label + 1
)
)
mask[coords] = False
log.debug(
"Update support for scale {:.2f}".format(self.kernels.scales[idx_scale])
)
data._support[idx_scale] |= mask
log.debug("Support sum: {}".format(data._support.sum()))
def _inverse_transform(self, data):
"""Do the inverse transform (reconstruct the image).
TODO: describe better what this does.
Parameters
----------
data : `~gammapy.detect.CWTData`
Images for inverse transform.
"""
data._transform_2d = np.sum(data._support * data._transform_3d, axis=0)
log.debug("Update model")
data._model += data._transform_2d * (data._transform_2d > 0)
log.debug("Model sum: {:.4f}".format(data._model.sum()))
log.debug("Model max: {:.4f}".format(data._model.max()))
def _is_converged(self, data):
"""Check if the algorithm has converged on current iteration.
TODO: document metric used, but not super important.
Parameters
----------
data : `~gammapy.detect.CWTData`
Images after iteration.
Returns
-------
answer : boolean
Answer if CWT has converged.
"""
log.debug("Check the convergence")
residual = data._counts - (data._model + data._approx)
variance = residual.var()
log.info("Residual sum: {:.4f}".format(residual.sum()))
log.info("Residual max: {:.4f}".format(residual.max()))
log.info("Residual variance: {:.4f}".format(residual.var()))
if self.previous_variance is None:
self.previous_variance = variance
return False
variance_ratio = abs(
(self.previous_variance - variance) / self.previous_variance
)
log.info("Variance ratio: {:.7f}".format(variance_ratio))
self.previous_variance = variance
return variance_ratio < self.tol
[docs] def analyze(self, data):
"""Run iterative filter peak algorithm.
The algorithm modifies the original data.
Parameters
----------
data : `~gammapy.detect.CWTData`
Input images.
"""
if self.history is not None:
self.history = [copy.deepcopy(data)]
for n_iter in range(self.max_iter):
log.info("************ Start iteration {} ************".format(n_iter + 1))
self._execute_iteration(data=data)
if self.history is not None:
log.debug("Save current data")
self.history.append(copy.deepcopy(data))
converge_answer = self._is_converged(data=data)
if converge_answer:
break
if converge_answer:
log.info("Convergence reached at iteration {}".format(n_iter + 1))
else:
log.info(
"Convergence not formally reached at iteration {}".format(n_iter + 1)
)
[docs]class CWTKernels:
"""Conduct arrays of kernels and scales for CWT algorithm.
Parameters
----------
n_scale : int
Number of scales.
min_scale : float
First scale used.
step_scale : float
Base scaling factor.
old : boolean (default False)
DEBUG attribute. If False, use astropy MaxicanHat kernels for kernel_base.
Attributes
----------
n_scale : int
Number of scales considered.
min_scale : float
First scale used.
step_scale : float
Base scaling factor.
scales : `~numpy.ndarray`
Grid of scales.
kern_base : dict
Dictionary of scale powers as keys and 2D kernel arrays.
(mexican hat) as values
kern_approx : `~numpy.ndarray`
2D Gaussian kernel array from maximum scale.
Examples
--------
>>> from gammapy.detect import CWTKernels
>>> kernels = CWTKernels(n_scale=3, min_scale=2.0, step_scale=2.6)
>>> print (kernels.info_table)
Name Source
---------------------------------- ----------------------
Number of scales 3
Minimal scale 2.0
Step scale 2.6
Scales [ 2. 5.2 13.52]
Kernels approx width 280
Kernels approx sum 0.99986288557
Kernels approx max 0.00012877518599
Kernels base width for 2.0 scale 40
Kernels base sum for 2.0 scale 0.000305108917065
Kernels base max for 2.0 scale 0.0315463182128
Kernels base width for 5.2 scale 107
Kernels base sum for 5.2 scale 0.000158044776015
Kernels base max for 5.2 scale 0.0050152112595
Kernels base width for 13.52 scale 280
Kernels base sum for 13.52 scale 0.000137114430344
Kernels base max for 13.52 scale 0.000740731187317
"""
def __init__(self, n_scale, min_scale, step_scale, old=False):
self.n_scale = n_scale
self.min_scale = min_scale
self.step_scale = step_scale
self.scales = np.array(
[min_scale * step_scale ** _ for _ in range(n_scale)], dtype=float
)
self.kern_base = {}
for idx_scale, scale in enumerate(self.scales):
if old:
self.kern_base[idx_scale] = difference_of_gauss_kernel(
scale, step_scale
)
else:
self.kern_base[idx_scale] = MexicanHat2DKernel(scale * step_scale).array
max_scale = min_scale * step_scale ** n_scale
self.kern_approx = Gaussian2DKernel(max_scale).array
def _info(self):
"""Return information about the object as a dict.
Returns
-------
info : dict
Information about object with str characteristic as keys and
characteristic results as values.
"""
info_dict = {}
info_dict["Number of scales"] = self.n_scale
info_dict["Minimal scale"] = self.min_scale
info_dict["Step scale"] = self.step_scale
info_dict["Scales"] = str(self.scales)
info_dict["Kernels approx width"] = len(self.kern_approx)
info_dict["Kernels approx sum"] = self.kern_approx.sum()
info_dict["Kernels approx max"] = self.kern_approx.max()
for idx_scale, scale in enumerate(self.scales):
info_dict["Kernels base width for {} scale".format(scale)] = len(
self.kern_base[idx_scale]
)
info_dict["Kernels base sum for {} scale".format(scale)] = self.kern_base[
idx_scale
].sum()
info_dict["Kernels base max for {} scale".format(scale)] = self.kern_base[
idx_scale
].max()
return info_dict
@property
def info_table(self):
"""Summary info table about the object.
Returns
-------
table : `~astropy.table.Table`
Information about the object.
"""
info_dict = self._info()
rows = []
for name in info_dict:
rows.append({"Name": name, "Source": info_dict[name]})
return Table(rows=rows, names=["Name", "Source"])
[docs]class CWTData:
"""Images for CWT algorithm.
Contains also input counts and background.
Parameters
----------
counts : `~gammapy.maps.WcsNDMap`
2D counts image.
background : `~gammapy.maps.WcsNDMap`
2D background image.
n_scale : int
Number of scales.
Examples
--------
>>> from gammapy.maps import Map
>>> from gammapy.detect import CWTData
>>> counts = Map.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz")
>>> background = Map.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-background.fits.gz")
>>> data = CWTData(counts=counts, background=background, n_scale=2)
"""
def __init__(self, counts, background, n_scale):
self._counts = np.array(counts.data, dtype=float)
self._background = np.array(background.data, dtype=float)
self._geom2d = counts.geom.copy()
scale_axis = MapAxis(np.arange(n_scale + 1))
self._geom3d = WcsGeom(
wcs=counts.geom.wcs, npix=counts.geom.npix, axes=[scale_axis]
)
shape_2d = self._counts.shape
self._model = np.zeros(shape_2d)
self._approx = np.zeros(shape_2d)
self._approx_bkg = np.zeros(shape_2d)
self._transform_2d = np.zeros(shape_2d)
shape_3d = n_scale, shape_2d[0], shape_2d[1]
self._transform_3d = np.zeros(shape_3d)
self._error = np.zeros(shape_3d)
self._support = np.zeros(shape_3d, dtype=bool)
@property
def counts(self):
"""2D counts input image (`~gammapy.maps.WcsNDMap`)."""
return WcsNDMap(self._geom2d, self._counts)
@property
def background(self):
"""2D background input image (`~gammapy.maps.WcsNDMap`)."""
return WcsNDMap(self._geom2d, self._background)
@property
def model(self):
"""2D model image (`~gammapy.maps.WcsNDMap`).
Positive version of transform_2d image.
Primordial initialized by zero array.
"""
return WcsNDMap(self._geom2d, self._model)
@property
def approx(self):
"""2D approx ??? image (`~gammapy.maps.WcsNDMap`).
In the course of iterations updated by convolution of
``counts - model - background`` with ``kern_approx``
Primordial initialized by zero array.
"""
return WcsNDMap(self._geom2d, self._approx)
@property
def approx_bkg(self):
"""2D approx bkg image (`~gammapy.maps.WcsNDMap`).
In the course of iterations updated by convolution of ``background`` with ``kern_approx``.
Primordial initialized by zero array.
"""
return WcsNDMap(self._geom2d, self._approx_bkg)
@property
def transform_2d(self):
"""2D transform ??? image (`~gammapy.maps.WcsNDMap`).
Created from transform_3d by summarize values per 0 axes.
Primordial initialized by zero array.
"""
return WcsNDMap(self._geom2d, self._transform_2d)
@property
def support_2d(self):
"""2D cube exclusion mask (`~gammapy.maps.WcsNDMap`).
Created from support_3d by OR-operation per 0 axis.
"""
support_2d = self._support.sum(0) > 0
return WcsNDMap(self._geom2d, support_2d)
@property
def residual(self):
"""2D residual image (`~gammapy.maps.WcsNDMap`).
Calculate as ``counts - model - approx``.
"""
residual = self._counts - (self._model + self._approx)
return WcsNDMap(self._geom2d, residual)
@property
def model_plus_approx(self):
"""TODO: document what this is."""
return WcsNDMap(self._geom2d, self._model + self._approx)
@property
def transform_3d(self):
"""3D transform ??? cube (`~gammapy.maps.WcsNDMap`).
Primordial initialized by zero array. In the course of
iterations updated by convolution of ``counts - total_background`` with kernel
for each scale (``total_background = model + background + approx``).
"""
return WcsNDMap(self._geom3d, self._transform_3d)
@property
def error(self):
"""3D error cube (`~gammapy.maps.WcsNDMap`).
Primordial initialized by zero array.
In the course of iterations updated by convolution of ``total_background``
with kernel^2 for each scale.
"""
return WcsNDMap(self._geom3d, self._error)
@property
def support_3d(self):
"""3D support (exclusion) cube (`~gammapy.maps.WcsNDMap`).
Primordial initialized by zero array.
"""
return WcsNDMap(self._geom3d, self._support)
@property
def max_scale_image(self):
"""Maximum scale image (`~gammapy.maps.WcsNDMap`)."""
# Previous version:
# idx_scale_max = np.argmax(self._transform_3d, axis=0)
# return kernels.scales[idx_scale_max] * (self._support.sum(0) > 0)
transform_2d_max = np.max(self._transform_3d, axis=0)
maximal_image = transform_2d_max * self.support_2d.data
return WcsNDMap(self._geom2d, maximal_image)
def __sub__(self, other):
data = CWTData(
counts=self.counts,
background=self.background,
n_scale=len(self._transform_3d),
)
data._model = self._model - other._model
data._approx = self._approx - other._approx
data._approx_bkg = self._approx_bkg - other._approx_bkg
data._transform_2d = self._transform_2d - other._transform_2d
data._transform_3d = self._transform_3d - other._transform_3d
data._error = self._error - other._error
data._support = self._support ^ other._support
return data
[docs] def images(self):
"""All the images in a dict.
Returns
-------
images : dict
Dictionary with keys {'counts', 'background', 'model', 'approx',
'approx_bkg', 'transform_2d', 'maximal', 'support_2d'}
and 2D `~numpy.ndarray` images as values.
"""
return dict(
counts=self.counts,
background=self.background,
model=self.model,
approx=self.approx,
approx_bkg=self.approx_bkg,
transform_2d=self.transform_2d,
model_plus_approx=self.model_plus_approx,
residual=self.residual,
maximal=self.max_scale_image,
support_2d=self.support_2d,
)
[docs] def cubes(self):
"""All the cubes in a dict.
Returns
-------
cubes : dict
Dictionary with keys {'transform_3d', 'error', 'support_3d'} and 3D
`~numpy.ndarray` cubes as values.
"""
return dict(
transform_3d=self.transform_3d, error=self.error, support=self.support_3d
)
@staticmethod
def _metrics_info(data, name):
"""Compute variance, mean, find max and min values and compute sum for given data.
Parameters
----------
data : `~numpy.ndarray`
2D image or 3D cube.
name : string
Name of the data.
Returns
-------
info : dict
The information about the data.
"""
return {
"Name": name,
"Shape": "2D image" if len(data.shape) == 2 else "3D cube",
"Variance": data.var(),
"Mean": data.mean(),
"Max value": data.max(),
"Min value": data.min(),
"Sum values": data.sum(),
}
[docs] def image_info(self, name):
"""Compute image info.
Compute variance, mean, find max and min values and compute sum for image with given name.
Return that information about the image.
Parameters
----------
name : string
Name of the image. Name can be as one of the follow: {'counts', 'background',
'model', 'approx', 'approx_bkg', 'transform_2d', 'model_plus_approx', 'residual',
'maximal', 'support_2d'}
Returns
-------
table : `~astropy.table.Table`
Information about the object.
"""
if name not in self.images():
raise ValueError(
"Incorrect name of image. It should be one of the following:"
"{'counts', 'background', 'model', 'approx', 'approx_bkg', "
"'transform_2d', 'model_plus_approx', 'residual', 'support_2d', "
"'maximal'}"
)
image = self.images()[name]
info_dict = self._metrics_info(data=image.data, name=name)
rows = []
for metric in info_dict:
rows.append({"Metrics": metric, "Source": info_dict[metric]})
return Table(rows=rows, names=["Metrics", "Source"])
[docs] def cube_info(self, name, per_scale=False):
"""Compute cube info.
Compute variance, mean, find max and min values and compute sum for image with given name.
Return that information about the image.
Parameters
----------
name : string
Name of the image. Name can be as one of the follow: {'transform_3d', 'error', 'support'}
per_scale : boolean, optional (default False)
If True, return information about the cube per all the scales.
Returns
-------
table : `~astropy.Table`
Information about the object.
"""
if name not in self.cubes():
raise ValueError(
"Incorrect name of cube. It should be one of the following:"
"{'transform_3d', 'error', 'support'}"
)
cube = self.cubes()[name]
rows = []
if per_scale:
mask = []
for index in range(len(cube.data)):
info_dict = self._metrics_info(data=cube.data[index], name=name)
for metric in info_dict:
rows.append(
{
"Scale power": index + 1,
"Metrics": metric,
"Source": info_dict[metric],
}
)
# For missing values in `Power scale` column
scale_mask = np.ones(len(info_dict), dtype=bool)
scale_mask[0] = False
mask.extend(scale_mask)
columns = ["Scale power", "Metrics", "Source"]
table = Table(rows=rows, names=columns, masked=True)
table["Scale power"].mask = mask
elif per_scale is False:
info_dict = self._metrics_info(data=cube.data, name=name)
for metric in info_dict:
rows.append({"Metrics": metric, "Source": info_dict[metric]})
columns = ["Metrics", "Source"]
table = Table(rows=rows, names=columns)
else:
raise ValueError("Incorrect value for per_scale attribute.")
return table
@property
def info_table(self):
"""Information about all the images and cubes.
Returns
-------
table : `~astropy.Table`
Information about the object.
"""
rows = []
for name, image in self.images().items():
info_dict = self._metrics_info(data=image.data, name=name)
rows.append(info_dict)
for name, cube in self.cubes().items():
info_dict = self._metrics_info(data=cube.data, name=name)
rows.append(info_dict)
columns = rows[0].keys()
return Table(rows=rows, names=columns)
[docs] def write(self, filename, overwrite=False):
"""Save results to FITS file.
Parameters
----------
filename : str
Fits file name.
overwrite : bool, optional (default False)
If True, overwrite file with name as filename.
"""
header = self._geom2d.make_header()
hdu_list = fits.HDUList()
hdu_list.append(fits.PrimaryHDU())
hdu_list.append(fits.ImageHDU(data=self._counts, header=header, name="counts"))
hdu_list.append(
fits.ImageHDU(data=self._background, header=header, name="background")
)
hdu_list.append(fits.ImageHDU(data=self._model, header=header, name="model"))
hdu_list.append(fits.ImageHDU(data=self._approx, header=header, name="approx"))
hdu_list.append(
fits.ImageHDU(data=self._transform_2d, header=header, name="transform_2d")
)
hdu_list.append(
fits.ImageHDU(data=self._approx_bkg, header=header, name="approx_bkg")
)
hdu_list.writeto(filename, overwrite=overwrite)