Source code for photutils.detection.core

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
This module implements the base class and star finder kernel for
detecting stars in an astronomical image.

Each star-finding class should define a method called ``find_stars``
that finds stars in an image.
"""

import abc
import math
import warnings

import numpy as np
from astropy.stats import gaussian_fwhm_to_sigma

from photutils.detection.peakfinder import find_peaks
from photutils.utils.exceptions import NoDetectionsWarning

__all__ = ['StarFinderBase']


[docs] class StarFinderBase(metaclass=abc.ABCMeta): """ Abstract base class for star finders. """
[docs] def __call__(self, data, mask=None): return self.find_stars(data, mask=mask)
@staticmethod def _find_stars(convolved_data, kernel, threshold, *, min_separation=0.0, mask=None, exclude_border=False): """ Find stars in an image. Parameters ---------- convolved_data : 2D array_like The convolved 2D array. kernel : `_StarFinderKernel` or 2D `~numpy.ndarray` The convolution kernel. ``StarFinder`` inputs the kernel as a 2D array. threshold : float The absolute image value above which to select sources. This threshold should be the threshold input to the star finder class multiplied by the kernel relerr. If ``convolved_data`` is a `~astropy.units.Quantity` array, then ``threshold`` must have the same units. min_separation : float, optional The minimum separation for detected objects in pixels. mask : 2D bool array, optional A boolean mask with the same shape as ``data``, where a `True` value indicates the corresponding element of ``data`` is masked. Masked pixels are ignored when searching for stars. exclude_border : bool, optional Set to `True` to exclude sources found within half the size of the convolution kernel from the image borders. The default is `False`, which is the mode used by IRAF's `DAOFIND <https://iraf.net/irafhelp.php?val=daofind>`_ and `starfind <https://iraf.net/irafhelp.php?val=starfind>`_. Returns ------- result : Nx2 `~numpy.ndarray` A Nx2 array containing the (x, y) pixel coordinates. """ # define a local footprint for the peak finder if min_separation == 0.0: # DAOStarFinder if isinstance(kernel, np.ndarray): footprint = np.ones(kernel.shape) else: footprint = kernel.mask.astype(bool) else: # define a local circular footprint for the peak finder idx = np.arange(-min_separation, min_separation + 1) xx, yy = np.meshgrid(idx, idx) footprint = np.array((xx**2 + yy**2) <= min_separation**2, dtype=int) # define the border exclusion region if exclude_border: if isinstance(kernel, np.ndarray): yborder = (kernel.shape[0] - 1) // 2 xborder = (kernel.shape[1] - 1) // 2 else: yborder = kernel.yradius xborder = kernel.xradius border_width = (yborder, xborder) else: border_width = None # find local peaks in the convolved data # suppress any NoDetectionsWarning from find_peaks with warnings.catch_warnings(): warnings.filterwarnings('ignore', category=NoDetectionsWarning) tbl = find_peaks(convolved_data, threshold, footprint=footprint, mask=mask, border_width=border_width) if tbl is None: return None return np.transpose((tbl['x_peak'], tbl['y_peak']))
[docs] @abc.abstractmethod def find_stars(self, data, mask=None): """ Find stars in an astronomical image. Parameters ---------- data : 2D array_like The 2D image array. mask : 2D bool array, optional A boolean mask with the same shape as ``data``, where a `True` value indicates the corresponding element of ``data`` is masked. Masked pixels are ignored when searching for stars. Returns ------- table : `~astropy.table.Table` or `None` A table of found stars. If no stars are found then `None` is returned. """ raise NotImplementedError('Needs to be implemented in a subclass.')
class _StarFinderKernel: """ Container class for a 2D Gaussian density enhancement kernel. The kernel has negative wings and sums to zero. It is used by both `DAOStarFinder` and `IRAFStarFinder`. Parameters ---------- fwhm : float The full-width half-maximum (FWHM) of the major axis of the Gaussian kernel in units of pixels. ratio : float, optional The ratio of the minor and major axis standard deviations of the Gaussian kernel. ``ratio`` must be strictly positive and less than or equal to 1.0. The default is 1.0 (i.e., a circular Gaussian kernel). theta : float, optional The position angle (in degrees) of the major axis of the Gaussian kernel, measured counter-clockwise from the positive x axis. sigma_radius : float, optional The truncation radius of the Gaussian kernel in units of sigma (standard deviation) [``1 sigma = FWHM / 2.0*sqrt(2.0*log(2.0))``]. The default is 1.5. normalize_zerosum : bool, optional Whether to normalize the Gaussian kernel to have zero sum, The default is `True`, which generates a density-enhancement kernel. Notes ----- The class attributes include the dimensions of the elliptical kernel and the coefficients of a 2D elliptical Gaussian function expressed as: ``f(x,y) = A * exp(-g(x,y))`` where ``g(x,y) = a*(x-x0)**2 + 2*b*(x-x0)*(y-y0) + c*(y-y0)**2`` References ---------- .. [1] https://en.wikipedia.org/wiki/Gaussian_function """ def __init__(self, fwhm, *, ratio=1.0, theta=0.0, sigma_radius=1.5, normalize_zerosum=True): if fwhm < 0: raise ValueError('fwhm must be positive.') if ratio <= 0 or ratio > 1: raise ValueError('ratio must be positive and less or equal ' 'than 1.') if sigma_radius <= 0: raise ValueError('sigma_radius must be positive.') self.fwhm = fwhm self.ratio = ratio self.theta = theta self.sigma_radius = sigma_radius self.xsigma = self.fwhm * gaussian_fwhm_to_sigma self.ysigma = self.xsigma * self.ratio theta_radians = np.deg2rad(self.theta) cost = np.cos(theta_radians) sint = np.sin(theta_radians) xsigma2 = self.xsigma**2 ysigma2 = self.ysigma**2 self.a = (cost**2 / (2.0 * xsigma2)) + (sint**2 / (2.0 * ysigma2)) # CCW self.b = 0.5 * cost * sint * ((1.0 / xsigma2) - (1.0 / ysigma2)) self.c = (sint**2 / (2.0 * xsigma2)) + (cost**2 / (2.0 * ysigma2)) # find the extent of an ellipse with radius = sigma_radius*sigma; # solve for the horizontal and vertical tangents of an ellipse # defined by g(x,y) = f self.f = self.sigma_radius**2 / 2.0 denom = (self.a * self.c) - self.b**2 # nx and ny are always odd # minimum kernel size is 5x5 self.nx = 2 * int(max(2, math.sqrt(self.c * self.f / denom))) + 1 self.ny = 2 * int(max(2, math.sqrt(self.a * self.f / denom))) + 1 self.xc = self.xradius = self.nx // 2 self.yc = self.yradius = self.ny // 2 # define the kernel on a 2D grid yy, xx = np.mgrid[0:self.ny, 0:self.nx] self.circular_radius = np.sqrt((xx - self.xc)**2 + (yy - self.yc)**2) self.elliptical_radius = (self.a * (xx - self.xc)**2 + 2.0 * self.b * (xx - self.xc) * (yy - self.yc) + self.c * (yy - self.yc)**2) self.mask = np.where( (self.elliptical_radius <= self.f) | (self.circular_radius <= 2.0), 1, 0).astype(int) self.npixels = self.mask.sum() # NOTE: the central (peak) pixel of gaussian_kernel has a value of 1.0 self.gaussian_kernel_unmasked = np.exp(-self.elliptical_radius) self.gaussian_kernel = self.gaussian_kernel_unmasked * self.mask # The denom represents (variance * npixels) denom = ((self.gaussian_kernel**2).sum() - (self.gaussian_kernel.sum()**2 / self.npixels)) self.relerr = 1.0 / np.sqrt(denom) # normalize the kernel to zero sum if normalize_zerosum: self.data = ((self.gaussian_kernel - (self.gaussian_kernel.sum() / self.npixels)) / denom) * self.mask else: # pragma: no cover self.data = self.gaussian_kernel self.shape = self.data.shape def _validate_brightest(brightest): """ Validate the ``brightest`` parameter. It must be >0 and an integer. """ if brightest is not None: if brightest <= 0: raise ValueError('brightest must be >= 0') bright_int = int(brightest) if bright_int != brightest: raise ValueError('brightest must be an integer') brightest = bright_int return brightest