# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
This module implements the IRAFStarFinder class.
"""
import inspect
import warnings
import astropy.units as u
import numpy as np
from astropy.nddata import extract_array
from astropy.table import QTable
from astropy.utils import lazyproperty
from photutils.detection.core import (StarFinderBase, _StarFinderKernel,
_validate_brightest)
from photutils.utils._convolution import _filter_data
from photutils.utils._misc import _get_meta
from photutils.utils._moments import _moments, _moments_central
from photutils.utils._quantity_helpers import isscalar, process_quantities
from photutils.utils.exceptions import NoDetectionsWarning
__all__ = ['IRAFStarFinder']
[docs]
class IRAFStarFinder(StarFinderBase):
"""
Detect stars in an image using IRAF's "starfind" algorithm.
`IRAFStarFinder` searches images for local density maxima that
have a peak amplitude greater than ``threshold`` above the local
background and have a PSF full-width at half-maximum similar to the
input ``fwhm``. The objects' centroid, roundness (ellipticity), and
sharpness are calculated using image moments.
Parameters
----------
threshold : float
The absolute image value above which to select sources.
If the star finder is run on an image that is a
`~astropy.units.Quantity` array, then ``threshold`` must have
the same units.
fwhm : float
The full-width half-maximum (FWHM) of the 2D circular Gaussian
kernel in units of pixels.
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))``].
minsep_fwhm : float, optional
The separation (in units of ``fwhm``) for detected objects. The
minimum separation is calculated as ``int((fwhm * minsep_fwhm) +
0.5)`` and is clipped to a minimum value of 2. Note that large
values may result in long run times.
sharplo : float, optional
The lower bound on sharpness for object detection.
sharphi : float, optional
The upper bound on sharpness for object detection.
roundlo : float, optional
The lower bound on roundness for object detection.
roundhi : float, optional
The upper bound on roundness for object detection.
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 starfind.
brightest : int, None, optional
The number of brightest objects to keep after sorting the source
list by flux. If ``brightest`` is set to `None`, all objects
will be selected.
peakmax : float, None, optional
The maximum allowed peak pixel value in an object. Only objects
whose peak pixel values are strictly smaller than ``peakmax``
will be selected. This may be used, for example, to exclude
saturated sources. If the star finder is run on an image that is
a `~astropy.units.Quantity` array, then ``peakmax`` must have
the same units. If ``peakmax`` is set to `None`, then no peak
pixel value filtering will be performed.
xycoords : `None` or Nx2 `~numpy.ndarray`, optional
The (x, y) pixel coordinates of the approximate centroid
positions of identified sources. If ``xycoords`` are input, the
algorithm will skip the source-finding step.
min_separation : `None` or float, optional
The minimum separation (in pixels) for detected objects. If
`None` then ``minsep_fwhm`` will be used, otherwise this keyword
overrides ``minsep_fwhm``. Note that large values may result in
long run times.
See Also
--------
DAOStarFinder
Notes
-----
If the star finder is run on an image that is a
`~astropy.units.Quantity` array, then ``threshold`` and ``peakmax``
must have the same units as the image.
For the convolution step, this routine sets pixels beyond the image
borders to 0.0. The equivalent parameters in IRAF's starfind are
``boundary='constant'`` and ``constant=0.0``.
IRAF's starfind uses ``hwhmpsf``, ``fradius``, and ``sepmin`` as
input parameters. The equivalent input values for `IRAFStarFinder`
are:
* ``fwhm = hwhmpsf * 2``
* ``sigma_radius = fradius * sqrt(2.0*log(2.0))``
* ``minsep_fwhm = 0.5 * sepmin``
The main differences between `~photutils.detection.DAOStarFinder`
and `~photutils.detection.IRAFStarFinder` are:
* `~photutils.detection.IRAFStarFinder` always uses a 2D
circular Gaussian kernel, while
`~photutils.detection.DAOStarFinder` can use an elliptical
Gaussian kernel.
* `IRAFStarFinder` internally calculates a "sky" background level
based on unmasked pixels within the kernel footprint.
* `~photutils.detection.IRAFStarFinder` calculates the objects'
centroid, roundness, and sharpness using image moments.
"""
def __init__(self, threshold, fwhm, sigma_radius=1.5, minsep_fwhm=2.5,
sharplo=0.5, sharphi=2.0, roundlo=0.0, roundhi=0.2,
exclude_border=False, brightest=None, peakmax=None,
xycoords=None, min_separation=None):
# here we validate the units, but do not strip them
inputs = (threshold, peakmax)
names = ('threshold', 'peakmax')
_ = process_quantities(inputs, names)
if not isscalar(threshold):
raise TypeError('threshold must be a scalar value.')
if not np.isscalar(fwhm):
raise TypeError('fwhm must be a scalar value.')
self.threshold = threshold
self.fwhm = fwhm
self.sigma_radius = sigma_radius
self.minsep_fwhm = minsep_fwhm
self.sharplo = sharplo
self.sharphi = sharphi
self.roundlo = roundlo
self.roundhi = roundhi
self.exclude_border = exclude_border
self.brightest = _validate_brightest(brightest)
self.peakmax = peakmax
if xycoords is not None:
xycoords = np.asarray(xycoords)
if xycoords.ndim != 2 or xycoords.shape[1] != 2:
raise ValueError('xycoords must be shaped as a Nx2 array')
self.xycoords = xycoords
self.kernel = _StarFinderKernel(self.fwhm, ratio=1.0, theta=0.0,
sigma_radius=self.sigma_radius)
if min_separation is not None:
if min_separation < 0:
raise ValueError('min_separation must be >= 0')
self.min_separation = min_separation
else:
self.min_separation = max(2, int((self.fwhm * self.minsep_fwhm)
+ 0.5))
def _get_raw_catalog(self, data, *, mask=None):
convolved_data = _filter_data(data, self.kernel.data, mode='constant',
fill_value=0.0,
check_normalization=False)
if self.xycoords is None:
xypos = self._find_stars(convolved_data, self.kernel,
self.threshold,
min_separation=self.min_separation,
mask=mask,
exclude_border=self.exclude_border)
else:
xypos = self.xycoords
if xypos is None:
warnings.warn('No sources were found.', NoDetectionsWarning)
return None
return _IRAFStarFinderCatalog(data, convolved_data, xypos, self.kernel,
sharplo=self.sharplo,
sharphi=self.sharphi,
roundlo=self.roundlo,
roundhi=self.roundhi,
brightest=self.brightest,
peakmax=self.peakmax)
[docs]
def find_stars(self, data, mask=None):
"""
Find stars in an astronomical image.
Parameters
----------
data : 2D array_like
The 2D image array. The image should be
background-subtracted.
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.QTable` or `None`
A table of found stars. `None` is returned if no stars are
found. The table contains the following parameters:
* ``id``: unique object identification number.
* ``xcentroid, ycentroid``: object centroid.
* ``fwhm``: object FWHM.
* ``sharpness``: object sharpness.
* ``roundness``: object roundness.
* ``pa``: object position angle (degrees counter clockwise from
the positive x axis).
* ``npix``: the total number of (positive) unmasked pixels.
* ``peak``: the peak, sky-subtracted, pixel value of the object.
* ``flux``: the object instrumental flux calculated as the
sum of sky-subtracted data values within the kernel
footprint.
* ``mag``: the object instrumental magnitude calculated as
``-2.5 * log10(flux)``.
"""
inputs = (data, self.threshold, self.peakmax)
names = ('data', 'threshold', 'peakmax')
_ = process_quantities(inputs, names)
cat = self._get_raw_catalog(data, mask=mask)
if cat is None:
return None
# apply all selection filters
cat = cat.apply_all_filters()
if cat is None:
return None
# create the output table
return cat.to_table()
class _IRAFStarFinderCatalog:
"""
Class to create a catalog of the properties of each detected star,
as defined by IRAF's ``starfind`` task.
Parameters
----------
data : 2D `~numpy.ndarray`
The 2D image. The image should be background-subtracted.
convolved_data : 2D `~numpy.ndarray`
The convolved 2D image. If ``data`` is a
`~astropy.units.Quantity` array, then ``convolved_data`` must
have the same units.
xypos : Nx2 `~numpy.ndarray`
A Nx2 array of (x, y) pixel coordinates denoting the central
positions of the stars.
kernel : `_StarFinderKernel`
The convolution kernel. This kernel must match the kernel used
to create the ``convolved_data``.
sharplo : float, optional
The lower bound on sharpness for object detection.
sharphi : float, optional
The upper bound on sharpness for object detection.
roundlo : float, optional
The lower bound on roundness for object detection.
roundhi : float, optional
The upper bound on roundness for object detection.
brightest : int, None, optional
The number of brightest objects to keep after sorting the source
list by flux. If ``brightest`` is set to `None`, all objects
will be selected.
peakmax : float, None, optional
The maximum allowed peak pixel value in an object. Only objects
whose peak pixel values are strictly smaller than ``peakmax``
will be selected. This may be used, for example, to exclude
saturated sources. If the star finder is run on an image that is
a `~astropy.units.Quantity` array, then ``peakmax`` must have
the same units. If ``peakmax`` is set to `None`, then no peak
pixel value filtering will be performed.
"""
def __init__(self, data, convolved_data, xypos, kernel, *, sharplo=0.2,
sharphi=1.0, roundlo=-1.0, roundhi=1.0, brightest=None,
peakmax=None):
# here we validate the units, but do not strip them
inputs = (data, convolved_data, peakmax)
names = ('data', 'convolved_data', 'peakmax')
_ = process_quantities(inputs, names)
self.data = data
unit = data.unit if isinstance(data, u.Quantity) else None
self.unit = unit
self.convolved_data = convolved_data
self.xypos = xypos
self.kernel = kernel
self.sharplo = sharplo
self.sharphi = sharphi
self.roundlo = roundlo
self.roundhi = roundhi
self.brightest = brightest
self.peakmax = peakmax
self.id = np.arange(len(self)) + 1
self.cutout_shape = kernel.shape
self.default_columns = ('id', 'xcentroid', 'ycentroid', 'fwhm',
'sharpness', 'roundness', 'pa', 'npix',
'peak', 'flux', 'mag')
def __len__(self):
return len(self.xypos)
def __getitem__(self, index):
# NOTE: we allow indexing/slicing of scalar (self.isscalar = True)
# instances in order to perform catalog filtering even for
# a single source
newcls = object.__new__(self.__class__)
# copy these attributes to the new instance
init_attr = ('data', 'unit', 'convolved_data', 'kernel',
'sharplo', 'sharphi', 'roundlo', 'roundhi', 'brightest',
'peakmax', 'cutout_shape', 'default_columns')
for attr in init_attr:
setattr(newcls, attr, getattr(self, attr))
# xypos determines ordering and isscalar
# NOTE: always keep as a 2D array, even for a single source
attr = 'xypos'
value = getattr(self, attr)[index]
setattr(newcls, attr, np.atleast_2d(value))
# index/slice the remaining attributes
keys = set(self.__dict__.keys()) & set(self._lazyproperties)
keys.add('id')
for key in keys:
value = self.__dict__[key]
# do not insert lazy attributes that are always scalar (e.g.,
# isscalar), i.e., not an array/list for each source
if np.isscalar(value):
continue
# value is always at least a 1D array, even for a single source
value = np.atleast_1d(value[index])
newcls.__dict__[key] = value
return newcls
@lazyproperty
def isscalar(self):
"""
Whether the instance is scalar (e.g., a single source).
"""
return self.xypos.shape == (1, 2)
@property
def _lazyproperties(self):
"""
Return all lazyproperties (even in superclasses).
"""
def islazyproperty(obj):
return isinstance(obj, lazyproperty)
return [i[0] for i in inspect.getmembers(self.__class__,
predicate=islazyproperty)]
def reset_ids(self):
"""
Reset the ID column to be consecutive integers.
"""
self.id = np.arange(len(self)) + 1
@lazyproperty
def sky(self):
"""
Calculate the sky background level.
The local sky level is roughly estimated using the IRAF starfind
calculation as the average value in the non-masked regions
within the kernel footprint.
"""
skymask = ~self.kernel.mask.astype(bool) # 1=sky, 0=obj
nsky = np.count_nonzero(skymask)
axis = (1, 2)
if nsky == 0.0: # pragma: no cover
sky = (np.max(self.cutout_data_nosub, axis=axis)
- np.max(self.cutout_convdata, axis=axis))
else:
sky = (np.sum(self.cutout_data_nosub * skymask, axis=axis)
/ nsky)
if self.unit is not None:
sky <<= self.unit
return sky
def make_cutouts(self, data):
cutouts = []
for xpos, ypos in self.xypos:
cutouts.append(extract_array(data, self.cutout_shape, (ypos, xpos),
fill_value=0.0))
value = np.array(cutouts)
if self.unit is not None:
value <<= self.unit
return value
@lazyproperty
def cutout_data_nosub(self):
return self.make_cutouts(self.data)
@lazyproperty
def cutout_data(self):
data = ((self.cutout_data_nosub - self.sky[:, np.newaxis, np.newaxis])
* self.kernel.mask)
# IRAF starfind discards negative pixels
data[data < 0] = 0.0
return data
@lazyproperty
def cutout_convdata(self): # pragma: no cover
return self.make_cutouts(self.convolved_data)
@lazyproperty
def npix(self):
return np.count_nonzero(self.cutout_data, axis=(1, 2))
@lazyproperty
def moments(self):
return np.array([_moments(arr, order=1) for arr in self.cutout_data])
@lazyproperty
def cutout_centroid(self):
moments = self.moments
# ignore divide-by-zero RuntimeWarning
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
ycentroid = moments[:, 1, 0] / moments[:, 0, 0]
xcentroid = moments[:, 0, 1] / moments[:, 0, 0]
return np.transpose((ycentroid, xcentroid))
@lazyproperty
def cutout_xcentroid(self):
return np.transpose(self.cutout_centroid)[1]
@lazyproperty
def cutout_ycentroid(self):
return np.transpose(self.cutout_centroid)[0]
@lazyproperty
def cutout_xorigin(self):
return np.transpose(self.xypos)[0] - self.kernel.xradius
@lazyproperty
def cutout_yorigin(self):
return np.transpose(self.xypos)[1] - self.kernel.yradius
@lazyproperty
def xcentroid(self):
return self.cutout_xcentroid + self.cutout_xorigin
@lazyproperty
def ycentroid(self):
return self.cutout_ycentroid + self.cutout_yorigin
@lazyproperty
def peak(self):
peaks = [np.max(arr) for arr in self.cutout_data]
return u.Quantity(peaks) if self.unit is not None else np.array(peaks)
@lazyproperty
def flux(self):
fluxes = [np.sum(arr) for arr in self.cutout_data]
if self.unit is not None:
fluxes = u.Quantity(fluxes)
else:
fluxes = np.array(fluxes)
return fluxes
@lazyproperty
def mag(self):
flux = self.flux
if isinstance(flux, u.Quantity):
flux = flux.value
return -2.5 * np.log10(flux)
@lazyproperty
def moments_central(self):
moments = np.array([_moments_central(arr, center=(xcen_, ycen_),
order=2)
for arr, xcen_, ycen_ in
zip(self.cutout_data, self.cutout_xcentroid,
self.cutout_ycentroid, strict=True)])
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
return moments / self.moments[:, 0, 0][:, np.newaxis, np.newaxis]
@lazyproperty
def mu_sum(self):
return self.moments_central[:, 0, 2] + self.moments_central[:, 2, 0]
@lazyproperty
def mu_diff(self):
return self.moments_central[:, 0, 2] - self.moments_central[:, 2, 0]
@lazyproperty
def fwhm(self):
return 2.0 * np.sqrt(np.log(2.0) * self.mu_sum)
@lazyproperty
def roundness(self):
# ignore divide-by-zero RuntimeWarning
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
return (np.sqrt(self.mu_diff**2
+ 4.0 * self.moments_central[:, 1, 1]**2)
/ self.mu_sum)
@lazyproperty
def sharpness(self):
return self.fwhm / self.kernel.fwhm
@lazyproperty
def pa(self):
pa = np.rad2deg(0.5 * np.arctan2(2.0 * self.moments_central[:, 1, 1],
self.mu_diff))
return np.where(pa < 0, pa + 180, pa)
def apply_filters(self):
"""
Filter the catalog.
"""
# remove all non-finite values - consider these non-detections
attrs = ('xcentroid', 'ycentroid', 'sharpness', 'roundness', 'pa',
'sky', 'peak', 'flux')
mask = np.count_nonzero(self.cutout_data, axis=(1, 2)) > 1
for attr in attrs:
mask &= np.isfinite(getattr(self, attr))
newcat = self[mask]
if len(newcat) == 0:
warnings.warn('No sources were found.', NoDetectionsWarning)
return None
# filter based on sharpness, roundness, and peakmax
mask = ((newcat.sharpness > newcat.sharplo)
& (newcat.sharpness < newcat.sharphi)
& (newcat.roundness > newcat.roundlo)
& (newcat.roundness < newcat.roundhi))
if newcat.peakmax is not None:
mask &= (newcat.peak < newcat.peakmax)
newcat = newcat[mask]
if len(newcat) == 0:
warnings.warn('Sources were found, but none pass the sharpness, '
'roundness, or peakmax criteria',
NoDetectionsWarning)
return None
return newcat
def select_brightest(self):
"""
Sort the catalog by the brightest fluxes and select the top
brightest sources.
"""
newcat = self
if self.brightest is not None:
idx = np.argsort(self.flux)[::-1][:self.brightest]
newcat = self[idx]
return newcat
def apply_all_filters(self):
"""
Apply all filters, select the brightest, and reset the source
IDs.
"""
cat = self.apply_filters()
if cat is None:
return None
cat = cat.select_brightest()
cat.reset_ids()
return cat
def to_table(self, columns=None):
table = QTable()
table.meta.update(_get_meta()) # keep table.meta type
if columns is None:
columns = self.default_columns
for column in columns:
table[column] = getattr(self, column)
return table