# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Tools for calculating properties of sources defined by an Aperture.
"""
import inspect
import warnings
from concurrent.futures import ThreadPoolExecutor
from copy import copy, deepcopy
from functools import cached_property
from typing import NamedTuple
import astropy.units as u
import numpy as np
from astropy.nddata import NDData
from astropy.stats import (SigmaClip, biweight_location, biweight_midvariance,
mad_std)
from photutils.aperture import Aperture, SkyAperture, region_to_aperture
from photutils.aperture._batch_photometry import (FLAG_COL_BBOX_CLIPPED,
FLAG_COL_MASKED,
FLAG_COL_N_PIXELS,
FLAG_COL_NONFINITE_DATA,
FLAG_COL_NONFINITE_ERROR,
FLAG_COL_SEG,
FLAG_COL_UNCORRECTED,
FLAG_COL_VALID,
batch_aperture_sums)
from photutils.aperture._batch_stats import (batch_aperture_gather,
batch_biweight, batch_gini,
batch_mad, batch_mean_var,
batch_moments, batch_order_stats,
batch_sigma_clip_center,
batch_sigma_clip_sum,
batch_sort_values)
from photutils.aperture._common import (SCALAR_COLLAPSE_TYPES,
batch_inputs_supported,
batch_mask_plane,
batch_segmentation_arrays,
collapse_scalar_value, unpack_nddata,
validate_array, validate_mask_method)
from photutils.aperture._segmentation import (make_segmentation_exclusion,
process_segmentation_inputs)
from photutils.aperture.core import (PixelAperture, _aperture_metadata,
_update_method_subpixels_docstring)
from photutils.aperture.flags import (APERTURE_FLAGS, _counts_to_flag_bits,
decode_aperture_flags)
from photutils.morphology import gini as gini_func
from photutils.utils._deprecation import (create_empty_deprecated_qtable,
deprecated, deprecated_getattr,
deprecated_positional_kwargs)
from photutils.utils._misc import _get_meta
from photutils.utils._moments import _image_moments
from photutils.utils._parameters import validate_table_columns
from photutils.utils._quantity_helpers import process_quantities
__all__ = ['ApertureStats']
# Scale factor that converts the median absolute deviation to a robust
# estimate of the standard deviation (1 / scipy.stats.norm.ppf(0.75)).
# This must match ``astropy.stats.mad_std`` and the value in
# ``photutils.aperture._batch_stats``.
_MAD_STD_SCALE = 1.482602218505602
# Remove in 4.0
_DEPRECATED_ATTRIBUTES: dict = {
'covar_sigx2': 'covariance_xx',
'covar_sigxy': 'covariance_xy',
'covar_sigy2': 'covariance_yy',
'cxx': 'ellipse_cxx',
'cxy': 'ellipse_cxy',
'cyy': 'ellipse_cyy',
'data_sumcutout': 'data_sum_cutout',
'error_sumcutout': 'error_sum_cutout',
'get_id': 'select_id',
'get_ids': 'select_ids',
'semimajor_sigma': 'semimajor_axis',
'semiminor_sigma': 'semiminor_axis',
'xcentroid': 'x_centroid',
'ycentroid': 'y_centroid',
}
class _BatchGather(NamedTuple):
"""
Container for the fast Cython batch-driver results shared by
`ApertureStats._fast_gather` and `ApertureStats._fast_sum`.
Each instance populates only the fields relevant to its footprint;
the remaining fields are `None`:
* center-value gather (``_fast_gather``): ``values``, ``local_x``,
``local_y``, ``starts``, ``counts``, ``overlap``, and
``flag_counts``, plus ``sorted_values`` (the packed
ascending-sorted surviving values) when sigma clipping is
applied
* ``sum_method`` gather (``_fast_sum``): ``sum_aper``, ``var_aper``,
``sum_area``, ``starts``, ``overlap``, and ``flag_counts``,
plus the packed member buffers (``sum_values``, ``sum_fracs``,
``sum_errsq``, and ``sum_counts``) while sigma clipping is applied
"""
values: np.ndarray = None
local_x: np.ndarray = None
local_y: np.ndarray = None
starts: np.ndarray = None
counts: np.ndarray = None
sum_aper: np.ndarray = None
var_aper: np.ndarray = None
sum_area: np.ndarray = None
overlap: np.ndarray = None
sum_values: np.ndarray = None
sum_fracs: np.ndarray = None
sum_errsq: np.ndarray = None
sum_counts: np.ndarray = None
flag_counts: np.ndarray = None
sorted_values: np.ndarray = None
# Public attributes that are never collapsed to a scalar for a scalar
# instance because they describe the whole object rather than a single
# per-source value (see ``ApertureStats.__getattribute__``). Any new
# public attribute that is a list, tuple, ndarray, or SkyCoord but is
# not a per-source value must be added here, otherwise a length-1 value
# will be silently collapsed to its first element for a scalar instance.
_SCALAR_EXCLUDE = frozenset({'default_columns', 'isscalar', 'labels',
'n_positions', 'properties',
'segmentation_image'})
[docs]
@_update_method_subpixels_docstring
class ApertureStats:
# numpydoc ignore: PR01,PR02,PR04,PR07
"""
Class to create a catalog of statistics for pixels within an
aperture.
Note that this class returns the statistics of the input
``data`` values within the aperture. It does not convert data
in surface brightness units to flux or counts. Conversion from
surface-brightness units should be performed before using this
class.
Parameters
----------
data : 2D `~numpy.ndarray`, `~astropy.units.Quantity`, \
`~astropy.nddata.NDData`
The 2D array from which to calculate the source properties.
For accurate source properties, ``data`` should be
background-subtracted. Non-finite ``data`` values (NaN and inf)
are automatically masked.
aperture : `~photutils.aperture.Aperture` or supported `~regions.Region`
The aperture or region to apply to the data. The aperture
or region object may contain more than one position. If the
input ``aperture`` is a `~photutils.aperture.SkyAperture` or
`~regions.SkyRegion` object, then a WCS must be input using
the ``wcs`` keyword. Region objects are converted to aperture
objects.
error : 2D `~numpy.ndarray` or `~astropy.units.Quantity`, optional
The total error array corresponding to the input ``data``
array. ``error`` is assumed to include *all* sources of
error, including the Poisson error of the sources (see
`~photutils.utils.calc_total_error`). ``error`` must have
the same shape as the input ``data``. If ``data`` is a
`~astropy.units.Quantity` array then ``error`` must be a
`~astropy.units.Quantity` array (and vice versa) with identical
units. Non-finite ``error`` values (NaN and +/- inf) are not
automatically masked, unless they are at the same position of
non-finite values in the input ``data`` array. Such pixels can
be masked using the ``mask`` keyword.
mask : 2D `~numpy.ndarray` (bool), optional
A boolean mask with the same shape as ``data`` where a `True`
value indicates the corresponding element of ``data`` is masked.
Masked data are excluded from all calculations. Non-finite
values (NaN and inf) in the input ``data`` are automatically
masked.
wcs : WCS object or `None`, optional
A world coordinate system (WCS) transformation that
supports the `astropy shared interface for WCS
<https://docs.astropy.org/en/stable/wcs/wcsapi.html>`_ (e.g.,
`astropy.wcs.WCS`, `gwcs.wcs.WCS`). ``wcs`` is required if
the input ``aperture`` is a `~photutils.aperture.SkyAperture`
or `~regions.SkyRegion` object. If `None`, then all sky-based
properties will be set to `None`.
sigma_clip : `None` or `astropy.stats.SigmaClip` instance, optional
A `~astropy.stats.SigmaClip` object that defines the sigma
clipping parameters. If `None` then no sigma clipping will
be performed.
sum_method : {'exact', 'center', 'subpixel'}, optional
The method used to determine the pixel weights (the fraction
of the pixel area covered by the aperture). This method
is used only for calculating the ``sum``, ``sum_error``,
``sum_aper_area``, ``data_sum_cutout``, and ``error_sum_cutout``
properties. All other properties use the "center" aperture mask
method. The following methods are available:
<method_bullets>
<subpixels_description>
ddof : int, optional
The delta degrees of freedom used when computing the ``var``
and ``std`` properties. The divisor used in the calculation
is ``N - ddof``, where ``N`` is the number of unmasked pixels
within the aperture. The default is ``ddof=0``, which gives the
population variance and standard deviation. Use ``ddof=1`` to
obtain the sample (unbiased) variance and standard deviation.
This keyword affects only the ``var`` and ``std`` properties.
All other properties are unaffected, including the ``mean_err``
and ``median_err`` standard errors, which are always computed
using the sample standard deviation. Apertures with ``N <=
ddof`` unmasked pixels have an undefined ``var`` and ``std`` and
are set to NaN.
local_bkg : float, `~numpy.ndarray`, `~astropy.units.Quantity`, or `None`
The per-pixel local background values to subtract from the data
before performing measurements. If input as an array, the order
of ``local_bkg`` values corresponds to the order of the input
``aperture`` positions. ``local_bkg`` must have the same length
as the input ``aperture`` or must be a scalar value, which
will be broadcast to all apertures. If `None`, then no local
background subtraction is performed. If the input ``data`` has
units, then ``local_bkg`` must be a `~astropy.units.Quantity`
with the same units.
<segmentation_descriptions>
n_threads : int, optional
The number of threads to use to gather the aperture pixel
values, compute the aperture sums, and reduce the per-source
statistics (e.g., the sorts, order statistics, robust
estimators, and moments). The default is 1 (no multithreading).
When ``n_threads`` > 1, the aperture positions are divided
into chunks and processed concurrently. The per-source
results are independent, so they are identical to the
single-threaded computation. The underlying kernels release
the Python global interpreter lock (GIL), so multithreading
can significantly speed up the statistics for many aperture
positions. Multithreading is used only on the fast batch code
path. Otherwise, the computation is serial.
Notes
-----
``data`` should be background-subtracted for accurate source
properties. In addition to global background subtraction, local
background subtraction can be performed using the ``local_bkg``
keyword values.
`~regions.Region` objects are converted to `Aperture` objects using
the :func:`region_to_aperture` function.
The returned statistics are measured for the pixels within the
input aperture at its input position. This class does not change
the position of the input aperture. This class returns the centroid
value of the pixels within the input aperture, but the input
aperture is not recentered at the measured centroid position
when making the measurements. If desired, you can create a new
`Aperture` object using the measured centroid and then re-run
`~photutils.aperture.ApertureStats`.
All properties other than the sum-related ones described below are
calculated using the "center" aperture-mask method, which assigns
aperture weights of either 0 or 1, so the ``data`` pixel values
are used directly and without weighting. This choice reflects a
fundamental limitation because, unlike the mean or variance, order
statistics (``min``, ``max``, ``median``) and robust estimators
(``mad_std``, ``biweight_location``, ``biweight_midvariance``) have
no standard, unambiguous definition for pixels with fractional
(partial) aperture weights. Accordingly, these quantities cannot be
rigorously computed from a weighted aperture footprint.
The input ``sum_method`` and ``subpixels`` keywords are
used to determine the aperture-mask method only for the
sum-related properties: ``sum``, ``sum_err``, ``sum_aper_area``,
``data_sum_cutout``, and ``error_sum_cutout`` (also listed in the
``SUM_FOOTPRINT_PROPERTIES`` class attribute). All other properties,
including ``mean``, ``median``, ``std``, and the morphological
properties, always use the "center" aperture-mask method regardless
of ``sum_method``. The default is ``sum_method='exact'``, which
produces exact aperture-weighted photometry.
The `flags` property reports quality conditions from both the
"center"-method footprint and the ``sum_method`` footprint (see the
`flags` docstring for the per-bit details).
The calculated statistics are always float64, regardless of the
input ``data`` dtype (`~astropy.units.Quantity` values with float64
dtype if the input ``data`` has units).
Examples
--------
>>> from photutils.datasets import make_4gaussians_image
>>> from photutils.aperture import CircularAperture, ApertureStats
>>> data = make_4gaussians_image()
>>> aper = CircularAperture((150, 25), 8)
>>> aperstats = ApertureStats(data, aper)
>>> print(aperstats.x_centroid)
149.99080259251238
>>> print(aperstats.y_centroid)
24.97484633000507
>>> print(aperstats.centroid)
[149.99080259 24.97484633]
>>> print(aperstats.mean, aperstats.median)
47.76300955780609 31.913789514433084
>>> print(aperstats.std)
39.193655383492974
>>> print(aperstats.sum)
9286.709206410273
>>> print(aperstats.sum_aper_area)
201.0619298297468 pix2
>>> # More than one aperture position
>>> aper2 = CircularAperture(((150, 25), (90, 60)), 10)
>>> aperstats2 = ApertureStats(data, aper2)
>>> print(aperstats2.x_centroid)
[149.98470724 89.97893946]
>>> print(aperstats2.sum)
[10177.62548482 36653.97704059]
"""
# Properties computed using the ``sum_method`` aperture-mask
# footprint (set by the ``sum_method``/``subpixels`` keywords).
# All other properties use the "center" aperture-mask method.
SUM_FOOTPRINT_PROPERTIES = ('sum', 'sum_err', 'sum_aper_area',
'data_sum_cutout', 'error_sum_cutout')
# Cached properties that are not per-source sliceable: the packed
# gather buffers and their reductions. ``__getitem__`` drops these
# from the sliced object, which recomputes them lazily from its
# sliced inputs. Any new cached property backed by the packed batch
# buffers must be added here, otherwise slicing will attempt to
# index the packed buffer per source and fail or corrupt it.
_NON_SLICEABLE_CACHES = frozenset({
'_batch_inputs', '_fast_gather', '_fast_sum', '_sorted_values',
'_order_stats', '_mean_var', '_mad', '_biweight', '_gini',
'_fast_cutouts_center'})
def __init__(self, data, aperture, *, error=None, mask=None, wcs=None,
sigma_clip=None, sum_method='exact', subpixels=5, ddof=0,
local_bkg=None, segmentation_image=None, labels=None,
mask_method='none', n_threads=1):
if isinstance(data, NDData):
data, error, mask, wcs = unpack_nddata(data, error, mask, wcs)
inputs = (data, error, local_bkg)
names = ('data', 'error', 'local_bkg')
inputs, unit = process_quantities(inputs, names)
(data, error, local_bkg) = inputs
self._data = validate_array(data, 'data')
self._data_unit = unit
self._validate_aperture(aperture)
aperture_meta = _aperture_metadata(aperture) # use input aperture
if isinstance(aperture, SkyAperture) and wcs is None:
msg = 'A wcs is required when using a SkyAperture'
raise ValueError(msg)
# Convert region to aperture if necessary
if not isinstance(aperture, Aperture):
aperture = region_to_aperture(aperture)
self.aperture = aperture
data_shape = self._data.shape
self._error = validate_array(error, 'error', shape=data_shape)
self._mask = validate_array(mask, 'mask', shape=data_shape)
self._wcs = wcs
if sigma_clip is not None and not isinstance(sigma_clip, SigmaClip):
msg = 'sigma_clip must be a SigmaClip instance'
raise TypeError(msg)
self.sigma_clip = sigma_clip
# Validate the mask-method keywords here so that an invalid
# value is reported at construction rather than at the first
# access of a measured property, far from its cause.
validate_mask_method(sum_method, subpixels,
method_name='sum_method')
self.sum_method = sum_method
self.subpixels = subpixels
if (isinstance(ddof, bool)
or not isinstance(ddof, (int, np.integer)) or ddof < 0):
msg = 'ddof must be a non-negative integer'
raise ValueError(msg)
self.ddof = ddof
if not isinstance(n_threads, (int, np.integer)) or n_threads < 1:
msg = 'n_threads must be a positive integer'
raise ValueError(msg)
self.n_threads = int(n_threads)
self._local_bkg = np.zeros(self.n_positions) # no local bkg
if local_bkg is not None:
local_bkg = np.atleast_1d(local_bkg)
if local_bkg.ndim != 1:
msg = 'local_bkg must be a 1D array'
raise ValueError(msg)
n_local_bkg = len(local_bkg)
if n_local_bkg not in (1, self.n_positions):
msg = ('local_bkg must be scalar or have the same length '
'as the input aperture')
raise ValueError(msg)
local_bkg = np.broadcast_to(local_bkg, self.n_positions)
if np.any(~np.isfinite(local_bkg)):
msg = ('local_bkg must not contain any non-finite '
'(e.g., inf or NaN) values')
raise ValueError(msg)
self._local_bkg = local_bkg # always an iterable
self._ids = np.arange(self.n_positions) + 1
self.default_columns = ['id', 'x_centroid', 'y_centroid',
'sky_centroid', 'sum', 'sum_err',
'sum_aper_area', 'center_aper_area',
'min', 'max', 'mean', 'median', 'mode',
'std', 'mad_std', 'var', 'biweight_location',
'biweight_midvariance', 'fwhm',
'semimajor_axis', 'semiminor_axis',
'orientation', 'eccentricity', 'flags']
self.meta = _get_meta()
self.meta.update(aperture_meta)
# Validate the segmentation-masking inputs and resolve the
# per-aperture source labels.
self.segmentation_image = segmentation_image
self.labels = labels
self.mask_method = mask_method
seg_positions = np.atleast_2d(self._pixel_aperture.positions)
(self._segmentation,
self._seg_labels) = process_segmentation_inputs(
segmentation_image, labels, mask_method,
seg_positions, self._data.shape)
@staticmethod
def _validate_aperture(aperture):
try:
from regions import Region
aper_types = (Aperture, Region)
except ImportError:
aper_types = Aperture
if not isinstance(aperture, aper_types):
msg = 'aperture must be an Aperture or Region object'
raise TypeError(msg)
return aperture
@property
def _cached_properties(self):
"""
A list of all class cached properties (even in superclasses).
The result is cached on the class to avoid repeated
introspection via `inspect.getmembers`.
"""
cls = self.__class__
attr = '_cached_properties_cache'
# Subclasses get their own cached-property list
if attr not in cls.__dict__:
def is_cached_property(obj):
return isinstance(obj, cached_property)
setattr(cls, attr,
[i[0] for i in inspect.getmembers(
cls, predicate=is_cached_property)])
return getattr(cls, attr)
@property
def properties(self):
"""
A sorted list of the built-in source properties.
"""
cached_properties = [name for name in self._cached_properties
if not name.startswith('_')]
# isscalar and n_positions are scalar values for the whole
# object, not per-source values, so they are not valid table
# columns
cached_properties.remove('isscalar')
cached_properties.remove('n_positions')
cached_properties.sort()
return cached_properties
def __getitem__(self, index):
if self.isscalar:
msg = (f'A scalar {self.__class__.__name__!r} object cannot '
'be indexed')
raise TypeError(msg)
newcls = object.__new__(self.__class__)
# Attributes defined in __init__ that are copied directly to the
# new class
init_attr = ('_data', '_data_unit', '_error', '_mask', '_wcs',
'sigma_clip', 'sum_method', 'subpixels', 'ddof',
'n_threads', 'default_columns', 'meta',
'_segmentation', 'segmentation_image', 'mask_method')
for attr in init_attr:
setattr(newcls, attr, getattr(self, attr))
# aperture determines isscalar (needed below)
newcls.aperture = self.aperture[index]
# Keep _ids as a 1D array so a scalar instance's id is always
# backed by a length-1 iterable (see the id property).
newcls._ids = np.atleast_1d(self._ids[index])
# Slice the per-aperture segmentation labels. Both the input
# ``labels`` and the resolved ``_seg_labels`` have one entry per
# aperture, so the sliced object reports the labels of the
# apertures it actually contains.
if self.labels is None:
newcls.labels = None
else:
newcls.labels = np.atleast_1d(self.labels)[index]
if self._seg_labels is None:
newcls._seg_labels = None
else:
newcls._seg_labels = np.atleast_1d(self._seg_labels[index])
# Slice evaluated cached-property objects
keys = set(self.__dict__.keys()) & set(self._cached_properties)
keys.add('_local_bkg') # iterable defined in __init__
# The packed gather buffers and their reductions are not
# per-source sliceable; the sliced object recomputes them
# lazily from its sliced inputs.
keys -= self._NON_SLICEABLE_CACHES
for key in keys:
value = self.__dict__[key]
# Do not insert attributes that are always scalar (e.g.,
# isscalar, n_positions), i.e., not an array/list for each
# source
if np.isscalar(value):
continue
try:
# Keep most _<attrs> as length-1 iterables
if (newcls.isscalar and key.startswith('_')
and key != '_pixel_aperture'):
if isinstance(value, np.ndarray):
val = value[:, np.newaxis][index]
else:
val = [value[index]]
else:
val = value[index]
except TypeError:
# Apply fancy indices (e.g., array/list or bool mask) to
# lists.
# See https://numpy.org/doc/stable/release/1.20.0-notes.html
# #arraylike-objects-which-do-not-define-len-and-getitem
arr = np.empty(len(value), dtype=object)
arr[:] = list(value)
val = arr[index].tolist()
newcls.__dict__[key] = val
return newcls
def __str__(self):
cls_name = f'<{self.__class__.__module__}.{self.__class__.__name__}>'
return f'{cls_name}\nLength: {self.n_positions}'
def __repr__(self):
return self.__str__()
def __len__(self):
if self.isscalar:
msg = f'Scalar {self.__class__.__name__!r} object has no len()'
raise TypeError(msg)
return self.n_positions
def __iter__(self):
for item in range(len(self)):
yield self.__getitem__(item)
# Remove in 4.0
def __getattr__(self, name):
return deprecated_getattr(self, name, _DEPRECATED_ATTRIBUTES,
since='3.0', until='4.0')
def __getattribute__(self, name):
# Collapse the leading position axis of the public per-source
# output attributes to a scalar when a single scalar aperture
# position is input (e.g., ``CircularAperture((10, 20), r=5)``).
# The scalar conversion is applied centrally here instead of
# being repeated on each individual property (see ``_array`` for
# the array-preserving accessor used by the table-building and
# flag-decoding machinery).
value = super().__getattribute__(name)
if (not name.startswith('_')
and name not in _SCALAR_EXCLUDE
and isinstance(value, SCALAR_COLLAPSE_TYPES)
and self.isscalar):
return collapse_scalar_value(value)
return value
def _array(self, name):
"""
Return the per-source attribute ``name`` with a leading (source)
axis, so that the table-building and flag-decoding machinery
always operates on arrays regardless of whether the instance is
scalar.
For a non-scalar instance this returns the attribute unchanged.
For a scalar instance the value is first read through the normal
attribute access, which applies the scalar collapse performed by
`__getattribute__` (e.g., ``moments`` becomes a ``(4, 4)`` array
and ``centroid`` a length-2 array). A single leading length-1
axis is then restored (``ndarray`` values are reshaped, while
other objects are wrapped in a length-1 list).
Reading the collapsed value rather than the raw stored value
is important. A per-source value can be stored either with
its leading length-1 axis intact (when computed directly
for a scalar instance) or with that axis already removed
(when a multi-source value is cached and then sliced to a
scalar via ``__getitem__``). Going through the scalar collapse
normalizes both cases to the same shape, so this method is
robust regardless of how the value was produced.
"""
value = getattr(self, name)
if not self.isscalar:
return value
if isinstance(value, np.ndarray):
return value[np.newaxis, ...]
return [value]
[docs]
@cached_property
def isscalar(self):
"""
Whether the instance is scalar (e.g., a single aperture
position).
"""
return self._pixel_aperture.isscalar
[docs]
def copy(self):
"""
Return a deep copy of this object.
Returns
-------
result : `ApertureStats`
A deep copy of this object.
"""
return deepcopy(self)
@cached_property
def _null_object(self):
"""
Return `None` values.
"""
return np.array([None] * self.n_positions)
@cached_property
def _null_value(self):
"""
Return np.nan values.
"""
values = np.empty(self.n_positions)
values.fill(np.nan)
return values
@property
def id(self):
"""
The aperture identification number(s).
"""
return self._ids
@property
@deprecated('3.1', alternative="the 'id' attribute", until='4.0')
def ids(self):
"""
The aperture identification number(s).
.. deprecated:: 3.1
Use the `id` attribute instead.
"""
return self.id
[docs]
def select_id(self, id_num):
"""
Return a new `ApertureStats` object for the input ID number
only.
Parameters
----------
id_num : int
The aperture ID number.
Returns
-------
result : `ApertureStats`
A new `ApertureStats` object containing only the source with
the input ID number.
"""
return self.select_ids(id_num)
[docs]
def select_ids(self, id_nums):
"""
Return a new `ApertureStats` object for the input ID numbers
only.
Parameters
----------
id_nums : list, tuple, or `~numpy.ndarray` of int
The aperture ID number(s).
Returns
-------
result : `ApertureStats`
A new `ApertureStats` object containing only the sources with
the input ID numbers.
Raises
------
TypeError
If this is a scalar `ApertureStats` object, which cannot be
indexed.
ValueError
If any input ID number is not a valid source ID number.
"""
if self.isscalar:
msg = (f'A scalar {self.__class__.__name__!r} object cannot '
'be indexed')
raise TypeError(msg)
for id_num in np.atleast_1d(id_nums):
if id_num not in self._array('id'):
msg = f'{id_num} is not a valid source ID number'
raise ValueError(msg)
sorter = np.argsort(self.id)
indices = sorter[np.searchsorted(self.id, id_nums, sorter=sorter)]
return self[indices]
[docs]
@deprecated_positional_kwargs(since='3.0', until='4.0')
def to_table(self, *, columns=None):
"""
Create a `~astropy.table.QTable` of source properties.
Parameters
----------
columns : str, list of str, `None`, optional
Names of columns, in order, to include in the output
`~astropy.table.QTable`. The allowed column names are any of
the `ApertureStats` properties. If ``columns`` is `None`,
then a default list of scalar-valued properties (as defined
by the ``default_columns`` attribute) will be used.
Returns
-------
table : `~astropy.table.QTable`
A table of sources properties with one row per source.
Raises
------
ValueError
If any name in ``columns`` is not a valid (or deprecated)
column name.
"""
if columns is None:
table_columns = self.default_columns
else:
# id is not included in self.properties because it is not
# a cached property
allowed_columns = set(self.properties) | set(self.default_columns)
# Remove 2D cutout images from the allowed columns
allowed_columns = {col for col in allowed_columns
if '_cutout' not in col}
deprecated_names = _DEPRECATED_ATTRIBUTES.copy()
# These are not valid column names, but are deprecated
# attributes that are still accessible in 3.x. They will be
# removed in 4.0.
invalid = ('data_sumcutout', 'error_sumcutout', 'get_id',
'get_ids')
for name in invalid:
deprecated_names.pop(name)
table_columns = validate_table_columns(
columns, allowed_columns, deprecated_names=deprecated_names)
# Replace with QTable in 4.0
tbl = create_empty_deprecated_qtable(
_DEPRECATED_ATTRIBUTES, since='3.0', until='4.0')
tbl.meta.update(self.meta) # keep tbl.meta type
for column in table_columns:
values = getattr(self, column)
# Column assignment requires an object with a length
if self.isscalar:
values = (values,)
# Use the canonical (non-deprecated) column name so that
# assigning into ``tbl`` does not trigger a second,
# redundant deprecation warning (the deprecated attribute
# access above already warned once).
canonical_column = _DEPRECATED_ATTRIBUTES.get(column, column)
tbl[canonical_column] = values
return tbl
[docs]
@cached_property
def n_positions(self):
"""
The number of positions for the input aperture.
"""
if self.isscalar:
return 1
return len(self._pixel_aperture)
@property
@deprecated('3.1', alternative="the 'n_positions' attribute",
until='4.0')
def n_apertures(self):
"""
The number of positions for the input aperture.
.. deprecated:: 3.1
Use the `n_positions` attribute instead.
"""
return self.n_positions
@cached_property
def _pixel_aperture(self):
"""
The input aperture as a PixelAperture.
"""
if isinstance(self.aperture, SkyAperture):
return self.aperture.to_pixel(self._wcs)
return self.aperture
@cached_property
def _batch_inputs(self):
"""
The validated inputs for the fast Cython batch driver, or `None`.
Returns a tuple of the contiguous arrays and scalar parameters
shared by the center-value and ``sum_method`` gathers (see
`_fast_gather` and `_fast_sum`), together with the fast
sigma-clip specification (`None` when no clipping is performed).
`None` is returned (and the mask-based code path is used) when
the aperture or inputs are not supported by the batch driver,
for example when an unsupported ``sigma_clip`` is used, the
aperture does not opt in to the batch driver, or the data/mask
dtypes are not supported. The set of supported inputs matches
`~photutils.aperture.PixelAperture._batch_photometry`.
"""
# A non-None sigma_clip is only supported when it maps to the
# fast clipping kernel; otherwise fall back to the mask path.
clip_spec = self._fast_clip_spec()
if self.sigma_clip is not None and clip_spec is None:
return None
aper = self._pixel_aperture
# Use the batch driver only if the aperture's own class
# opted in via the _enable_batch_photometry decorator (see
# PixelAperture._batch_photometry).
if type(aper)._batch_photometry_class is not type(aper):
return None
spec = aper._batch_shape_params()
if spec is None:
return None
data = self._data
error = self._error
if not batch_inputs_supported(data, error, self._mask):
return None
# Non-finite ``data`` values are always folded into the mask
# plane so the batch kernels can skip the per-pixel finiteness
# test (and defer the pixel-value load to contributing pixels
# only). This matches the mask-based path, which masks
# non-finite data before any segmentation correction.
mask = batch_mask_plane(data, self._mask, mask_nonfinite=True)
seg_arr, labels_arr, seg_code = batch_segmentation_arrays(
self._segmentation, self._seg_labels, self.mask_method)
shape_code, params = spec
sum_use_exact, sum_subpixels = aper._translate_mask_method(
self.sum_method, self.subpixels)
ext_x, ext_y = aper._xy_extents
off_x, off_y = aper._xy_bbox_offset
if error is not None:
error = np.ascontiguousarray(error, dtype=np.float64)
return (np.ascontiguousarray(data, dtype=np.float64), error, mask,
np.ascontiguousarray(aper._positions, dtype=np.float64),
shape_code, np.array(params, dtype=np.float64),
float(ext_x), float(ext_y), float(off_x), float(off_y),
sum_use_exact, sum_subpixels,
np.ascontiguousarray(self._local_bkg, dtype=np.float64),
seg_arr, labels_arr, seg_code, clip_spec)
@cached_property
def _fast_gather(self):
"""
The fast Cython "center"-method value gather, or `None`.
When supported, this walks each aperture bounding box once and
returns a packed buffer of the unmasked "center"-method pixel
values (and their cutout coordinates). The higher-tier
statistics (order statistics and image moments) are reduced
lazily from this buffer. The ``sum_method`` aperture sum, error
variance, and area are computed separately and lazily by
`_fast_sum`, so they are not paid for unless requested.
When a supported ``sigma_clip`` is set, the packed center buffer
is sigma-clipped per source before being returned, so the
downstream reductions operate on the clipped data transparently.
`None` is returned (and the mask-based code path is used) when
the fast batch driver is unavailable (see `_batch_inputs`).
The result is a `_BatchGather` with the center-value fields
(``values``, ``local_x``, ``local_y``, ``starts``, ``counts``),
``overlap``, and ``flag_counts`` populated. The ``sum_method``
fields are `None`.
When ``n_threads`` > 1, the positions are divided into chunks
that are gathered (and sigma clipped) concurrently and then
merged (see `_merge_center_gathers`).
"""
inputs = self._batch_inputs
if inputs is None:
return None
(data, _error, mask, positions, shape_code, params, ext_x, ext_y,
off_x, off_y, _sum_use_exact, _sum_subpixels, local_bkg, seg_arr,
labels_arr, seg_code, clip_spec) = inputs
def gather_chunk(pos, bkg, labels):
(values, lx, ly, starts, counts, overlap,
flag_counts) = batch_aperture_gather(
data, mask, pos, shape_code, params, ext_x, ext_y,
off_x, off_y, bkg, seg_arr, labels, seg_code)
gather = _BatchGather(values=values, local_x=lx, local_y=ly,
starts=starts, counts=counts,
overlap=overlap, flag_counts=flag_counts)
if clip_spec is not None:
gather = self._apply_center_clip(gather, clip_spec)
return gather
chunks = self._batch_chunks(positions, local_bkg, labels_arr)
if chunks is None:
return gather_chunk(positions, local_bkg, labels_arr)
with ThreadPoolExecutor(max_workers=len(chunks[0])) as executor:
gathers = list(executor.map(gather_chunk, *chunks))
return self._merge_center_gathers(gathers)
@cached_property
def _fast_sum(self):
"""
The fast Cython ``sum_method`` aperture gather, or `None`.
When supported, this walks each aperture bounding box once and
returns the ``sum_method`` aperture sum, error variance, and
area. This is computed only when the ``sum``, ``sum_err``, or
``sum_aper_area`` properties are requested, so the more
expensive exact/subpixel overlap fractions are not evaluated for
the value-statistics-only use cases.
When a supported ``sigma_clip`` is set, the packed ``sum_method``
members are sigma-clipped per source and the aperture sum,
variance, and area are recomputed over the survivors.
`None` is returned (and the mask-based code path is used) when
the fast batch driver is unavailable (see `_batch_inputs`).
The result is a `_BatchGather` with the ``sum_aper``,
``var_aper``, ``sum_area``, ``starts``, ``overlap``, and
``flag_counts`` fields populated. The center-value fields are
`None`.
When ``n_threads`` > 1, the positions are divided into chunks
that are computed (and sigma clipped) concurrently and then
merged. The merged result keeps only the flat per-source outputs
(see `_merge_sum_gathers`).
"""
inputs = self._batch_inputs
if inputs is None:
return None
(data, error, mask, positions, shape_code, params, ext_x, ext_y,
off_x, off_y, sum_use_exact, sum_subpixels, local_bkg, seg_arr,
labels_arr, seg_code, clip_spec) = inputs
emit_sum = 1 if clip_spec is not None else 0
def sum_chunk(pos, bkg, labels):
(sums, sum_var, area, overlap, starts, sum_values, sum_fracs,
sum_errsq, scounts, flag_counts) = batch_aperture_sums(
data, error, mask, pos, shape_code, params, ext_x, ext_y,
off_x, off_y, sum_use_exact, sum_subpixels, seg_arr,
labels, seg_code, bkg, emit_sum)
gather = _BatchGather(starts=starts, sum_aper=sums,
var_aper=sum_var, sum_area=area,
overlap=overlap, sum_values=sum_values,
sum_fracs=sum_fracs, sum_errsq=sum_errsq,
sum_counts=scounts,
flag_counts=flag_counts)
if clip_spec is not None:
gather = self._apply_sum_clip(gather, clip_spec)
return gather
chunks = self._batch_chunks(positions, local_bkg, labels_arr)
if chunks is None:
return sum_chunk(positions, local_bkg, labels_arr)
with ThreadPoolExecutor(max_workers=len(chunks[0])) as executor:
gathers = list(executor.map(sum_chunk, *chunks))
return self._merge_sum_gathers(gathers)
def _batch_chunks(self, positions, local_bkg, labels_arr):
"""
Split the per-source batch-driver inputs into per-thread chunks,
or return `None` for a single-chunk (serial) computation.
The number of chunks is ``min(n_threads, n_sources)``. Row
slices of the C-contiguous input arrays are themselves
C-contiguous, so the chunks can be passed directly to the Cython
drivers.
Returns
-------
chunks : tuple of list or `None`
A ``(positions, local_bkg, labels)`` tuple of per-chunk
lists, or `None` when only one chunk would be used.
"""
n_chunks = min(self.n_threads, positions.shape[0])
if n_chunks <= 1: # 0 for empty positions
return None
pos_chunks = np.array_split(positions, n_chunks)
bkg_chunks = np.array_split(local_bkg, n_chunks)
if labels_arr is None:
labels_chunks = [None] * n_chunks
else:
labels_chunks = np.array_split(labels_arr, n_chunks)
return pos_chunks, bkg_chunks, labels_chunks
@staticmethod
def _merge_center_gathers(gathers):
"""
Merge per-chunk center-value gathers into a single
`_BatchGather` equivalent to a single-chunk gather.
The packed buffers (``values``, ``local_x``, ``local_y``, and
``sorted_values`` when sigma clipping is applied) and the
per-source arrays are concatenated. The per-source ``starts``
are offset by the cumulative packed-buffer length of the
preceding chunks.
"""
starts = []
offset = 0
for gather in gathers:
starts.append(gather.starts + offset)
offset += gather.values.shape[0]
sorted_values = None
if gathers[0].sorted_values is not None:
sorted_values = np.concatenate(
[gather.sorted_values for gather in gathers])
return _BatchGather(
values=np.concatenate([g.values for g in gathers]),
local_x=np.concatenate([g.local_x for g in gathers]),
local_y=np.concatenate([g.local_y for g in gathers]),
starts=np.concatenate(starts),
counts=np.concatenate([g.counts for g in gathers]),
overlap=np.concatenate([g.overlap for g in gathers]),
flag_counts=np.concatenate([g.flag_counts for g in gathers]),
sorted_values=sorted_values)
@staticmethod
def _merge_sum_gathers(gathers):
"""
Merge per-chunk ``sum_method`` gathers into a single
`_BatchGather`.
Only the flat per-source outputs are kept. The packed member
buffers and their ``starts`` are dropped (`None`). They are
consumed within each chunk (by the sigma clipping) and no
downstream consumer reads them from the merged result.
"""
return _BatchGather(
sum_aper=np.concatenate([g.sum_aper for g in gathers]),
var_aper=np.concatenate([g.var_aper for g in gathers]),
sum_area=np.concatenate([g.sum_area for g in gathers]),
overlap=np.concatenate([g.overlap for g in gathers]),
flag_counts=np.concatenate([g.flag_counts for g in gathers]))
def _threaded_reduction(self, func, buffers, starts, counts,
per_source=()):
"""
Run a per-source packed-buffer reduction, dividing the sources
into chunks that are processed concurrently.
The kernel is called as ``func(*buffers, starts, counts,
*per_source)``. When ``n_threads`` > 1, the sources are
divided into contiguous ranges; each worker receives the
corresponding packed-buffer regions (with the ``starts``
rebased to the region) and per-source array slices, so the
kernels run concurrently on disjoint data. The per-chunk
outputs are concatenated, making the result identical to the
single-chunk call. A kernel returning a packed buffer (e.g.,
``batch_sort_values``) also concatenates correctly because the
chunk regions partition the buffer in order.
Parameters
----------
func : callable
The reduction kernel.
buffers : tuple of `~numpy.ndarray`
The packed per-pixel buffers (equal lengths), passed as
the leading kernel arguments.
starts, counts : `~numpy.ndarray`
The per-source start offsets into the buffers and the
per-source pixel counts.
per_source : tuple of `~numpy.ndarray`, optional
Additional per-source arrays, passed as the trailing
kernel arguments.
Returns
-------
result : `~numpy.ndarray` or tuple of `~numpy.ndarray`
The kernel output(s), concatenated over the chunks.
"""
n_src = len(counts)
n_chunks = min(self.n_threads, n_src)
if n_chunks <= 1: # 0 for empty positions
return func(*buffers, starts, counts, *per_source)
buffer_len = len(buffers[0])
src_edges = np.arange(n_chunks + 1) * n_src // n_chunks
def reduce_chunk(i0, i1):
buf0 = starts[i0]
buf1 = starts[i1] if i1 < n_src else buffer_len
chunk_buffers = [buffer[buf0:buf1] for buffer in buffers]
return func(*chunk_buffers, starts[i0:i1] - buf0,
counts[i0:i1],
*[arr[i0:i1] for arr in per_source])
with ThreadPoolExecutor(max_workers=n_chunks) as executor:
results = list(executor.map(reduce_chunk, src_edges[:-1],
src_edges[1:]))
if isinstance(results[0], tuple):
return tuple(np.concatenate(parts)
for parts in zip(*results, strict=True))
return np.concatenate(results)
def _fast_clip_spec(self):
"""
The fast sigma-clip parameters, or `None`.
Returns ``(sigma_lower, sigma_upper, maxiters, cenfunc_code,
stdfunc_code)`` when ``sigma_clip`` is a `SigmaClip` instance
supported by the fast clipping kernel, and `None` otherwise (in
which case the mask-based path is used). The fast path supports
the string ``cenfunc`` values 'median'/'mean'/'biweight', the
string ``stdfunc`` values 'std'/'mad_std'/'biweight', and no
spatial growing.
"""
sc = self.sigma_clip
if sc is None or not isinstance(sc, SigmaClip):
return None
if not (isinstance(sc.cenfunc, str)
and sc.cenfunc in ('median', 'mean', 'biweight')):
return None
if not (isinstance(sc.stdfunc, str)
and sc.stdfunc in ('std', 'mad_std', 'biweight')):
return None
if sc.grow: # False or 0 is supported; spatial growing is not
return None
maxiters = sc.maxiters
if maxiters is None or maxiters == np.inf:
maxiters = -1
else:
maxiters = int(maxiters)
cenfunc_code = {'median': 0, 'mean': 1, 'biweight': 2}[sc.cenfunc]
stdfunc_code = {'std': 0, 'mad_std': 1, 'biweight': 2}[sc.stdfunc]
return (float(sc.sigma_lower), float(sc.sigma_upper), maxiters,
cenfunc_code, stdfunc_code)
def _apply_center_clip(self, gather, clip_spec):
"""
Sigma-clip the packed center buffer and return a `_BatchGather`
whose center buffer reflects the per-source clipped data.
"""
sigma_lower, sigma_upper, maxiters, cenfunc, stdfunc = clip_spec
(cvalues, clx, cly, cstarts, ccounts,
csorted) = batch_sigma_clip_center(
gather.values, gather.local_x, gather.local_y, gather.starts,
gather.counts, sigma_lower, sigma_upper, maxiters, cenfunc,
stdfunc)
return gather._replace(values=cvalues, local_x=clx, local_y=cly,
starts=cstarts, counts=ccounts,
sorted_values=csorted)
def _apply_sum_clip(self, gather, clip_spec):
"""
Sigma-clip the packed ``sum_method`` member buffers and return
a `_BatchGather` whose aperture sum, variance, and area reflect
the per-source clipped data. The packed member buffers are
dropped (set to `None`) so their memory can be reclaimed.
"""
sigma_lower, sigma_upper, maxiters, cenfunc, stdfunc = clip_spec
has_error = 1 if self._error is not None else 0
(csum, cvar, carea) = batch_sigma_clip_sum(
gather.sum_values, gather.sum_fracs, gather.sum_errsq,
gather.starts, gather.sum_counts, sigma_lower, sigma_upper,
maxiters, cenfunc, stdfunc, has_error)
return gather._replace(sum_aper=csum, var_aper=cvar, sum_area=carea,
sum_values=None, sum_fracs=None,
sum_errsq=None, sum_counts=None)
@cached_property
def _sorted_values(self):
"""
The packed per-source ascending-sorted center pixel values, or
`None`.
When the fast gather path is active (see `_fast_gather`), each
source's packed pixel values are sorted once. The sorted buffer
is cached and shared by the order statistics (``min``, ``max``,
``median``), ``mad_std``, and the biweight estimators, so the
per-source sort is performed only once. When sigma clipping is
applied, the sorted surviving values produced by the clipping
kernel are reused directly (the clipping already sorts each
source's values to compute the clip bounds). `None` is returned
when the fast path is unavailable.
When ``n_threads`` > 1, the per-source sorts run concurrently
(see `_threaded_reduction`).
"""
gather = self._fast_gather
if gather is None:
return None
values, starts, counts = gather.values, gather.starts, gather.counts
if gather.sorted_values is not None:
return (gather.sorted_values, starts, counts)
sorted_values = self._threaded_reduction(
batch_sort_values, (values,), starts, counts)
return (sorted_values, starts, counts)
@cached_property
def _order_stats(self):
"""
The per-source ``(min, max, median)`` arrays, or `None`.
Reduced from the cached sorted buffer (see `_sorted_values`).
`None` when the fast path is unavailable.
"""
sorted_values = self._sorted_values
if sorted_values is None:
return None
values, starts, counts = sorted_values
return self._threaded_reduction(batch_order_stats, (values,),
starts, counts)
@cached_property
def _mean_var(self):
"""
The per-source ``(mean, var)`` arrays, or `None`.
The population (``ddof=0``) variance is returned. Computed
directly from the packed center buffer (no sort required).
`None` when the fast path is unavailable.
"""
gather = self._fast_gather
if gather is None:
return None
return self._threaded_reduction(batch_mean_var, (gather.values,),
gather.starts, gather.counts)
@cached_property
def _mad(self):
"""
The per-source unscaled median absolute deviation, or `None`.
Reduced from the cached sorted buffer (see `_sorted_values`).
`None` when the fast path is unavailable.
"""
sorted_values = self._sorted_values
if sorted_values is None:
return None
values, starts, counts = sorted_values
return self._threaded_reduction(batch_mad, (values,), starts, counts)
@cached_property
def _biweight(self):
"""
The per-source ``(biweight_location, biweight_midvariance)``
arrays, or `None`.
Reused from the cached sorted buffer, median (`_order_stats`),
and unscaled MAD (`_mad`). `None` when the fast path is
unavailable.
"""
sorted_values = self._sorted_values
if sorted_values is None:
return None
_, _, median = self._order_stats
values, starts, counts = sorted_values
return self._threaded_reduction(batch_biweight, (values,),
starts, counts,
per_source=(median, self._mad))
@cached_property
def _gini(self):
"""
The per-source Gini coefficient, or `None`.
Computed from the packed center buffer (the absolute values are
sorted internally). `None` when the fast path is unavailable.
"""
gather = self._fast_gather
if gather is None:
return None
return self._threaded_reduction(batch_gini, (gather.values,),
gather.starts, gather.counts)
def _finalize_value_stat(self, fast_result, stat_func, *,
square_unit=False, apply_unit=True):
"""
Return a per-source value statistic, using the fast Cython
reduction when available and otherwise the mask-based per-source
code path.
Parameters
----------
fast_result : `~numpy.ndarray` or `None`
The per-source statistic from the fast reduction, or `None`
to use the mask-based fallback.
stat_func : callable
The fallback callable applied to each source's 1D array of
unmasked pixel values when the fast path is unavailable.
square_unit : bool, optional
Whether the statistic has squared data units (e.g. variance).
apply_unit : bool, optional
Whether to apply the data unit at all (`False` for
dimensionless statistics such as the Gini coefficient).
"""
if fast_result is not None:
result = fast_result.copy()
else:
result = np.array([stat_func(arr)
for arr in self._data_values_center])
if apply_unit:
unit = self._data_unit
if unit is not None:
if square_unit:
unit = unit**2
result <<= unit
return result
@cached_property
def _aperture_masks_center(self):
"""
The aperture masks (`ApertureMask`) generated with the 'center'
method, always as an iterable.
"""
aperture_masks = self._pixel_aperture.to_mask(method='center')
if self.isscalar:
aperture_masks = (aperture_masks,)
return aperture_masks
@cached_property
def _aperture_masks(self):
"""
The aperture masks (`ApertureMask`) generated with the
``sum_method`` method, always as an iterable.
"""
aperture_masks = self._pixel_aperture.to_mask(method=self.sum_method,
subpixels=self.subpixels)
if self.isscalar:
aperture_masks = (aperture_masks,)
return aperture_masks
@cached_property
def _overlap_slices(self):
"""
The aperture mask overlap slices with the data, always as an
iterable.
The overlap slices are the same for all aperture mask methods.
"""
overlap_slices = []
for apermask in self._aperture_masks_center:
(slc_large, slc_small) = apermask.get_overlap_slices(
self._data.shape)
overlap_slices.append((slc_large, slc_small))
return overlap_slices
@cached_property
def _data_cutouts(self):
"""
The local-background-subtracted unmasked data cutouts using the
aperture bounding box, always as an iterable.
"""
cutouts = []
for (slices, local_bkg) in zip(self._overlap_slices,
self._local_bkg, strict=True):
if slices[0] is None:
cutout = None # no aperture overlap with the data
else:
# Copy is needed to preserve input data because masks are
# applied to these cutouts later
cutout = (self._data[slices[0]].astype(float, copy=True)
- local_bkg)
cutouts.append(cutout)
return cutouts
def _make_aperture_cutouts(self, aperture_masks, *, count_clipped=True):
"""
Make aperture-weighted cutouts for the data and variance, and
cutouts for the total mask and aperture mask weights.
Parameters
----------
aperture_masks : list of `ApertureMask`
A list of `ApertureMask` objects.
count_clipped : bool, optional
Whether to count the number of sigma-clipped pixels per
source. Only `_footprint_flag_inputs` for the "center"
footprint uses this count, so it is skipped (left at 0)
for the ``sum_method`` footprint to avoid the wasted
computation.
Returns
-------
result : list of `~numpy.ndarray`
A list of cutout arrays for the data, variance,
mask and weight arrays for each source (aperture
position), followed by the overlap indicator, per-source
flag counts (see the ``FLAG_COL_*`` constants in
`photutils.aperture._batch_photometry`, with semantics
identical to the batch drivers), and the number of
sigma-clipped pixels (always 0 when ``count_clipped`` is
`False`).
"""
# Use a local copy of the SigmaClip instance because SigmaClip
# stores internal state on the instance during calls. This
# method is reachable from two different cached properties (the
# center- and sum-footprint cutouts), which do not share a lock,
# so calling a shared instance from multiple threads could
# silently corrupt the results.
sigma_clip = copy(self.sigma_clip)
data_cutouts = []
variance_cutouts = []
mask_cutouts = []
weight_cutouts = []
overlaps = []
flag_counts = []
n_clipped = []
positions = np.atleast_2d(self._pixel_aperture.positions)
for idx, (data_cutout, apermask, slices) in enumerate(
zip(self._data_cutouts, aperture_masks,
self._overlap_slices, strict=True)):
fc_row = np.zeros(8, dtype=np.intp)
n_clipped_ = 0
slc_large, slc_small = slices
if slc_large is None: # aperture does not overlap the data
overlap = False
data_cutout = np.array([np.nan])
variance_cutout = np.array([np.nan])
mask_cutout = np.array([False])
weight_cutout = np.array([np.nan])
else:
# Create a mask of non-finite ``data`` values combined
# with the input ``mask`` array
nonfinite_mask = ~np.isfinite(data_cutout)
if self._mask is not None:
user_mask = self._mask[slc_large]
data_mask = nonfinite_mask | user_mask
else:
user_mask = None
data_mask = nonfinite_mask
error_cutout = (None if self._error is None
else self._error[slc_large])
# Apply segmentation-based masking and/or symmetric
# neighbor correction
exclude = None
affected = None
if (self._segmentation is not None
and self.mask_method != 'none'):
segm_cutout = self._segmentation[slc_large]
cutout_xycen = (positions[idx, 0] - slc_large[1].start,
positions[idx, 1] - slc_large[0].start)
(data_cutout, error_cutout, exclude,
affected) = make_segmentation_exclusion(
self.mask_method, segm_cutout,
self._seg_labels[idx], data=data_cutout,
error=error_cutout, base_mask=data_mask,
cutout_xycen=cutout_xycen)
pre_seg_mask = data_mask
data_mask = data_mask | exclude
overlap = True
aperweight_cutout = apermask.data[slc_small]
weight_cutout = aperweight_cutout * ~data_mask
# Per-source pixel counts for the quality flags,
# matching the batch-driver semantics (computed before
# any sigma clipping)
weighted = aperweight_cutout > 0
fc_row[FLAG_COL_N_PIXELS] = np.count_nonzero(weighted)
fc_row[FLAG_COL_BBOX_CLIPPED] = (
aperweight_cutout.shape != apermask.data.shape)
if user_mask is not None:
fc_row[FLAG_COL_MASKED] = np.count_nonzero(
weighted & user_mask)
fc_row[FLAG_COL_NONFINITE_DATA] = np.count_nonzero(
weighted & nonfinite_mask & ~user_mask)
else:
fc_row[FLAG_COL_NONFINITE_DATA] = np.count_nonzero(
weighted & nonfinite_mask)
if affected is not None:
pre_valid = weighted & ~pre_seg_mask
fc_row[FLAG_COL_SEG] = np.count_nonzero(
affected & pre_valid)
if self.mask_method == 'correct':
# In 'correct' mode, the excluded pixels are
# exactly the uncorrectable neighbor pixels
fc_row[FLAG_COL_UNCORRECTED] = np.count_nonzero(
exclude & pre_valid)
fc_row[FLAG_COL_VALID] = np.count_nonzero(
weighted & ~data_mask)
if error_cutout is not None:
fc_row[FLAG_COL_NONFINITE_ERROR] = np.any(
~np.isfinite(error_cutout) & weighted & ~data_mask)
# Apply the aperture mask; for "exact" and "subpixel"
# this is an expanded boolean mask using the aperture
# mask zero values
mask_cutout = (aperweight_cutout == 0) | data_mask
data_cutout = data_cutout.copy()
if sigma_clip is None:
# data_cutout will have zeros where mask_cutout is True
data_cutout *= ~mask_cutout
else:
# To input a mask, SigmaClip needs a MaskedArray
data_cutout_ma = np.ma.masked_array(data_cutout,
mask=mask_cutout)
data_sigclip = sigma_clip(data_cutout_ma)
# Define a mask of only the sigma-clipped pixels
sigclip_mask = data_sigclip.mask & ~mask_cutout
if count_clipped:
n_clipped_ = np.count_nonzero(sigclip_mask)
weight_cutout *= ~sigclip_mask
mask_cutout = data_sigclip.mask
data_cutout = data_sigclip.filled(0.0)
# Need to apply the aperture weights
data_cutout *= aperweight_cutout
if self._error is None:
variance_cutout = None
else:
# Apply the exact weights and total mask;
# error_cutout will have zeros where mask_cutout is True
variance = error_cutout**2
variance_cutout = (variance * aperweight_cutout**2
* ~mask_cutout)
data_cutouts.append(data_cutout)
variance_cutouts.append(variance_cutout)
mask_cutouts.append(mask_cutout)
weight_cutouts.append(weight_cutout)
overlaps.append(overlap)
flag_counts.append(fc_row)
n_clipped.append(n_clipped_)
# Use zip (instead of np.transpose) because these may contain
# arrays that have different shapes
return list(zip(data_cutouts, variance_cutouts, mask_cutouts,
weight_cutouts, overlaps, flag_counts, n_clipped,
strict=True))
@cached_property
def _aperture_cutouts_center(self):
"""
Aperture-weighted cutouts for the data, variance, total mask,
and aperture weights using the "center" aperture mask method.
"""
return self._make_aperture_cutouts(self._aperture_masks_center)
@cached_property
def _aperture_cutouts(self):
"""
Aperture-weighted cutouts for the data, variance, total mask,
and aperture weights using the input ``sum_method`` aperture
mask method.
"""
# The sigma-clipped pixel count is not tracked for this
# footprint: only the "center" footprint's flags use it (see
# `_footprint_flag_inputs`), so counting it here would be wasted
# computation.
return self._make_aperture_cutouts(self._aperture_masks,
count_clipped=False)
@cached_property
def _mask_cutout_center(self):
"""
Boolean mask cutouts representing the total mask.
The total mask is a combination of the input ``mask``,
non-finite ``data`` values, the cutout aperture mask using the
"center" method, and the sigma-clip mask.
"""
fast = self._fast_cutouts_center
if fast is not None:
return fast[2]
return list(zip(*self._aperture_cutouts_center, strict=True))[2]
@cached_property
def _mask_cutout(self):
"""
Boolean mask cutouts representing the total mask.
The total mask is a combination of the input ``mask``,
non-finite ``data`` values, the cutout aperture mask using the
``sum_method`` method, and the sigma-clip mask.
"""
return list(zip(*self._aperture_cutouts, strict=True))[2]
@cached_property
def _fast_cutouts_center(self):
"""
The center-method cutout arrays reconstructed from the packed
gather buffers, or `None`.
Returns a ``(data, variance, mask, weight)`` tuple of
per-source cutout lists identical to the mask-based
`_aperture_cutouts_center` values. The packed buffer holds
exactly the surviving pixels (unmasked, finite, inside the
center-method aperture, segmentation-kept, and sigma-clip
surviving) together with their cutout coordinates, so the total
mask is `True` for every pixel absent from the buffer, and the
data, variance, and weight cutouts are zero there (matching
the mask-based arrays). Sources whose bounding box does not
overlap the data get the same single-NaN sentinel arrays as the
mask-based path.
`None` is returned (and the mask-based path is used) when
the fast gather is unavailable or when segmentation neighbor
correction is applied: ``mask_method='correct'`` replaces
neighbor-pixel data *and error* values with their mirrored
values, and the mirrored error values are not recoverable from
the packed buffers.
"""
gather = self._fast_gather
if gather is None or (self._segmentation is not None
and self.mask_method == 'correct'):
return None
# The kernel's per-source cutout origin is the bounding box
# clipped to the data, using the same integer bounds as
# ``PixelAperture._bbox_bounds`` (upper bounds exclusive).
bounds = self._pixel_aperture._bbox_bounds
ny, nx = self._data.shape
x0 = np.maximum(bounds[:, 0], 0)
x1 = np.minimum(bounds[:, 1], nx)
y0 = np.maximum(bounds[:, 2], 0)
y1 = np.minimum(bounds[:, 3], ny)
overlaps = ((bounds[:, 0] < nx) & (bounds[:, 1] > 0)
& (bounds[:, 2] < ny) & (bounds[:, 3] > 0))
error = self._error
starts, counts = gather.starts, gather.counts
values = gather.values
local_x, local_y = gather.local_x, gather.local_y
data_cutouts = []
variance_cutouts = []
mask_cutouts = []
weight_cutouts = []
for idx, overlap in enumerate(overlaps):
if not overlap:
# Match the mask-based no-overlap sentinels
data_cutouts.append(np.array([np.nan]))
variance_cutouts.append(np.array([np.nan]))
mask_cutouts.append(np.array([False]))
weight_cutouts.append(np.array([np.nan]))
continue
shape = (y1[idx] - y0[idx], x1[idx] - x0[idx])
slc = (slice(y0[idx], y1[idx]), slice(x0[idx], x1[idx]))
i0 = starts[idx]
i1 = i0 + counts[idx]
lx = local_x[i0:i1]
ly = local_y[i0:i1]
data_cutout = np.zeros(shape)
data_cutout[ly, lx] = values[i0:i1]
if self.sigma_clip is None:
# The mask-based path multiplies the data by the
# boolean masks, so non-finite pixels are NaN in the
# cutout data instead of zero. (With sigma clipping
# it instead fills every masked pixel with zero.)
data_cutout[~np.isfinite(self._data[slc])] = np.nan
mask_cutout = np.ones(shape, dtype=bool)
mask_cutout[ly, lx] = False
weight_cutout = np.zeros(shape)
weight_cutout[ly, lx] = 1.0
if error is None:
variance_cutout = None
else:
# Multiplying the full variance cutout by the weights
# replicates the mask-based arithmetic exactly,
# including NaN for masked non-finite error values. The
# center-method weights are 0 or 1, so squaring them is
# not needed.
variance_cutout = error[slc] ** 2 * weight_cutout
data_cutouts.append(data_cutout)
variance_cutouts.append(variance_cutout)
mask_cutouts.append(mask_cutout)
weight_cutouts.append(weight_cutout)
return (data_cutouts, variance_cutouts, mask_cutouts,
weight_cutouts)
def _make_masked_array_center(self, array):
"""
Return a list of cutout masked arrays using the ``_mask_cutout``
mask.
Units are not applied.
"""
return [np.ma.masked_array(arr, mask=mask)
for arr, mask in zip(array, self._mask_cutout_center,
strict=True)]
def _make_masked_array(self, array):
"""
Return a list of cutout masked arrays using the
``_mask_cutout`` mask.
Units are not applied.
"""
return [np.ma.masked_array(arr, mask=mask)
for arr, mask in zip(array, self._mask_cutout, strict=True)]
[docs]
@cached_property
def data_cutout(self):
"""
A 2D aperture-weighted cutout from the data using the aperture
mask with the "center" method as a `~numpy.ma.MaskedArray`.
The cutout does not have units due to current limitations of
masked quantity arrays.
The mask is `True` for pixels from the input ``mask``,
non-finite ``data`` values (NaN and inf), sigma-clipped pixels
within the aperture, and pixels where the aperture mask has zero
weight.
"""
fast = self._fast_cutouts_center
if fast is not None:
return self._make_masked_array_center(fast[0])
return self._make_masked_array_center(
list(zip(*self._aperture_cutouts_center, strict=True))[0])
[docs]
@cached_property
def data_sum_cutout(self):
"""
A 2D aperture-weighted cutout from the data using the aperture
mask with the input ``sum_method`` method as a
`~numpy.ma.MaskedArray`.
The cutout does not have units due to current limitations of
masked quantity arrays.
The mask is `True` for pixels from the input ``mask``,
non-finite ``data`` values (NaN and inf), sigma-clipped pixels
within the aperture, and pixels where the aperture mask has zero
weight.
"""
return self._make_masked_array(list(zip(*self._aperture_cutouts,
strict=True))[0])
@cached_property
def _variance_cutout_center(self):
"""
A 2D variance cutout weighted by the squared aperture mask
weights, using the aperture mask with the input "center" method,
as a `~numpy.ma.MaskedArray`.
The cutout does not have units due to current limitations of
masked quantity arrays.
The mask is `True` for pixels from the input ``mask``,
non-finite ``data`` values (NaN and inf), sigma-clipped pixels
within the aperture, and pixels where the aperture mask has zero
weight.
"""
if self._error is None:
return self._null_object
fast = self._fast_cutouts_center
if fast is not None:
return self._make_masked_array_center(fast[1])
return self._make_masked_array_center(
list(zip(*self._aperture_cutouts_center, strict=True))[1])
@cached_property
def _variance_cutout(self):
"""
A 2D variance cutout weighted by the squared aperture mask
weights, using the aperture mask with the input ``sum_method``
method, as a `~numpy.ma.MaskedArray`.
The cutout does not have units due to current limitations of
masked quantity arrays.
The mask is `True` for pixels from the input ``mask``,
non-finite ``data`` values (NaN and inf), sigma-clipped pixels
within the aperture, and pixels where the aperture mask has zero
weight.
"""
if self._error is None:
return self._null_object
return self._make_masked_array(list(zip(*self._aperture_cutouts,
strict=True))[1])
[docs]
@cached_property
def error_sum_cutout(self):
"""
A 2D aperture-weighted error cutout using the aperture mask with
the input ``sum_method`` method as a `~numpy.ma.MaskedArray`.
The cutout values are the pixel errors multiplied by the
aperture mask weights, so the quadrature sum of the unmasked
cutout values equals `sum_err`.
The cutout does not have units due to current limitations of
masked quantity arrays.
The mask is `True` for pixels from the input ``mask``,
non-finite ``data`` values (NaN and inf), sigma-clipped pixels
within the aperture, and pixels where the aperture mask has zero
weight.
"""
if self._error is None:
return self._null_object
return [np.sqrt(var) for var in self._variance_cutout]
@cached_property
def _weight_cutout_center(self):
"""
A 2D `~numpy.ma.MaskedArray` cutout from the aperture mask
weights array using the aperture bounding box.
The aperture mask weights are for the "center" method.
The mask is `True` for pixels outside the aperture mask, pixels
from the input ``mask``, non-finite ``data`` values (NaN and
inf), and sigma-clipped pixels.
"""
fast = self._fast_cutouts_center
if fast is not None:
return self._make_masked_array_center(fast[3])
return self._make_masked_array_center(
list(zip(*self._aperture_cutouts_center, strict=True))[3])
@cached_property
def _weight_cutout(self):
"""
A 2D `~numpy.ma.MaskedArray` cutout from the aperture mask
weights array using the aperture bounding box.
The aperture mask weights are for the ``sum_method`` method.
The mask is `True` for pixels outside the aperture mask, pixels
from the input ``mask``, non-finite ``data`` values (NaN and
inf), and sigma-clipped pixels.
"""
return self._make_masked_array(list(zip(*self._aperture_cutouts,
strict=True))[3])
@cached_property
def _moment_data_cutout(self):
"""
A list of 2D `~numpy.ndarray` cutouts from the data.
Masked pixels are set to zero in these arrays (zeros do not
contribute to the image moments). The aperture mask weights are
for the "center" method.
These arrays are used to derive moment-based properties.
"""
data = deepcopy(self._array('data_cutout'))
cutouts = []
for arr in data:
if arr.size == 1 and np.isnan(arr[0]): # no aperture overlap
arr_ = np.empty((2, 2))
arr_.fill(np.nan)
else:
arr_ = arr.data
arr_[arr.mask] = 0.0
cutouts.append(arr_)
return cutouts
@cached_property
def _all_masked(self):
"""
True if all pixels within the aperture are masked.
"""
return np.array([np.all(mask) for mask in self._mask_cutout_center])
@cached_property
def _overlap(self):
"""
True if there is no overlap of the aperture with the data.
"""
return list(zip(*self._aperture_cutouts_center, strict=True))[4]
def _footprint_flag_inputs(self, footprint):
"""
The per-source flag inputs for the given aperture footprint.
Parameters
----------
footprint : {'center', 'sum'}
The aperture footprint: the "center" mask method used by the
value statistics, or the ``sum_method`` mask method used by
the sum properties.
Returns
-------
flag_counts : `~numpy.ndarray`
The per-source flag counts (see the ``FLAG_COL_*``
constants in `photutils.aperture._batch_photometry`).
overlap : `~numpy.ndarray` (bool)
Whether the aperture bounding box overlaps the data.
n_kept : `~numpy.ndarray` or `None`
For the center footprint, the per-source number of valid
pixels remaining after sigma clipping; `None` for the sum
footprint.
"""
if footprint == 'center':
gather = self._fast_gather
if gather is not None:
return (gather.flag_counts, np.asarray(gather.overlap),
np.asarray(gather.counts))
cutouts = self._aperture_cutouts_center
else:
gather = self._fast_sum
if gather is not None:
return gather.flag_counts, np.asarray(gather.overlap), None
cutouts = self._aperture_cutouts
flag_counts = np.array([cut[5] for cut in cutouts])
overlap = np.array([cut[4] for cut in cutouts])
n_kept = None
if footprint == 'center':
n_clipped = np.array([cut[6] for cut in cutouts])
n_kept = flag_counts[:, FLAG_COL_VALID] - n_clipped
return flag_counts, overlap, n_kept
def _footprint_flags(self, footprint):
"""
The per-source flag bits for the given aperture footprint.
The bounding-box-clipped candidate sources are resolved to the
precise outside-weight test using per-source aperture masks.
"""
flag_counts, overlap, _ = self._footprint_flag_inputs(footprint)
if footprint == 'center':
method, subpixels = 'center', 1
else:
method, subpixels = self.sum_method, self.subpixels
candidates = flag_counts[:, FLAG_COL_BBOX_CLIPPED].astype(bool)
w_out = self._pixel_aperture._resolve_outside_weights(
self._data.shape, method=method, subpixels=subpixels,
candidates=candidates)
return _counts_to_flag_bits(flag_counts, overlap, w_out)
@cached_property
def _base_flags(self):
"""
The count-based value-statistics flag bits (1D int array).
These are the "center"-method footprint bits plus the sigma-clip
and ``ddof`` bits. The `flags` property combines them with the
``sum_method`` footprint bits and the ``undefined_shape`` and
``singular_covariance`` bits.
"""
# The gather kernel and the center-method cutouts do not
# evaluate error values, so the non-finite-error bit is stripped
# here. `flags` picks it up from the sum footprint, the only
# footprint where error values are evaluated.
flags = (self._footprint_flags('center')
& ~APERTURE_FLAGS.NON_FINITE_ERROR)
flag_counts, _, n_kept = self._footprint_flag_inputs('center')
n_valid = flag_counts[:, FLAG_COL_VALID]
if self.sigma_clip is not None:
n_clipped = n_valid - n_kept
flags[n_clipped > 0] |= APERTURE_FLAGS.SIGMA_CLIPPED
flags[(n_valid > 0)
& (n_kept == 0)] |= APERTURE_FLAGS.ALL_CLIPPED
if self.ddof > 0:
w_in = flag_counts[:, FLAG_COL_N_PIXELS]
flags[(w_in > 0)
& (n_kept <= self.ddof)] |= APERTURE_FLAGS.TOO_FEW_PIXELS
return flags
@cached_property
def _undefined_shape_mask(self):
"""
Boolean mask (1D) marking sources whose net flux is not
positive.
The net flux is the zeroth image moment of the unmasked
"center"-method pixels. When it is zero or negative, the
centroid and the covariance-derived shape properties are
undefined or unreliable. Sources with no valid pixels (no
overlap, fully masked, or fully sigma clipped) are not flagged
here. They are already reported by the overlap, masking, and
clipping bits.
"""
m00 = self._array('moments')[:, 0, 0]
# NaN where a source has no valid pixels
n_pixels = self._center_n_pixels
return np.isfinite(m00) & (m00 <= 0) & np.isfinite(n_pixels)
@cached_property
def _singular_covariance_mask(self):
"""
Boolean mask (1D) marking sources with a singular or nearly
singular covariance matrix.
A source is flagged when the minor-axis variance (the smaller
eigenvalue of its normalized second-moment covariance matrix)
falls below ``1/12``. This flags both unresolved, nearly
point-like sources, where the covariance determinant drops below
``(1/12)**2``, and rank-1 degenerate sources, where one axis is
unresolved while the other is extended (which the determinant
alone would miss). Sources with undefined moments (no overlap or
fully masked) have a non-finite determinant and are not flagged
here. They are already reported by the overlap and masking bits.
"""
covar = self._raw_covariance
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
covar_det = np.linalg.det(covar)
# Smaller eigenvalue (the minor-axis variance) of each 2x2
# symmetric covariance matrix, via the closed form
# lambda = tr/2 -/+ sqrt((tr/2)**2 - det). The discriminant
# ((lambda1 - lambda2)/2)**2 is non-negative for a real
# symmetric matrix; clip tiny negative rounding to zero.
half_trace = 0.5 * (covar[:, 0, 0] + covar[:, 1, 1])
disc = np.maximum(half_trace**2 - covar_det, 0.0)
min_eigval = half_trace - np.sqrt(disc)
# ``1/12`` is the variance of a uniform distribution across a
# single pixel, and hence the smallest second moment a resolved
# source can have given finite pixel size. The determinant test
# (det < (1/12)**2) flags the isotropic case where both axes are
# unresolved; the eigenvalue test (min_eigval < 1/12)
# additionally flags rank-1 degeneracy (one unresolved axis) and
# covariance matrices that are not positive semidefinite
# (det < 0).
delta = 1.0 / 12
finite = np.isfinite(covar_det) & np.isfinite(min_eigval)
return finite & ((covar_det < delta**2) | (min_eigval < delta))
[docs]
@cached_property
@_update_method_subpixels_docstring
def flags(self):
# numpydoc ignore: RT01
"""
The bitwise quality flags.
The footprint-based flags (e.g., ``'masked_pixels'``,
``'non_finite_data'``, and the overlap flags) are evaluated
on the union of the "center"-method footprint used by the
value statistics and the ``sum_method`` footprint used by the
sum properties. The ``'non_finite_error'`` flag is evaluated
on the ``sum_method`` footprint. The ``'sigma_clipped'``,
``'all_clipped'``, and ``'too_few_pixels'`` flags are evaluated
on the value-statistics footprint. The ``'undefined_shape'`` and
``'singular_covariance'`` flags are always evaluated. Accessing
``flags`` computes the moment and covariance properties if they
have not already been computed (the results are cached and
shared with the corresponding shape properties).
See `~photutils.aperture.decode_aperture_flags` for decoding
flag values. The flags are:
<flag_descriptions>
"""
# The | allocates a new array, so the in-place |= below never
# mutates the cached `_base_flags`
flags = self._base_flags | self._footprint_flags('sum')
flags[self._undefined_shape_mask] |= (
APERTURE_FLAGS.UNDEFINED_SHAPE)
flags[self._singular_covariance_mask] |= (
APERTURE_FLAGS.SINGULAR_COVARIANCE)
return flags
[docs]
def decode_flags(self, *, return_bit_values=False):
"""
Decode the source quality flags into individual components.
This is a convenience method that calls
`~photutils.aperture.decode_aperture_flags` with the `flags`
property.
Parameters
----------
return_bit_values : bool, optional
If `True`, return the decoded bit flags (integers) instead
of the flag names (strings).
Returns
-------
decoded : dict
A dictionary mapping each aperture position `id` to the list
of its active flag names (or bit values). The entries follow
the position order.
See Also
--------
photutils.aperture.decode_aperture_flags
Examples
--------
>>> import numpy as np
>>> from photutils.aperture import ApertureStats, CircularAperture
>>> data = np.ones((25, 25))
>>> mask = np.zeros(data.shape, dtype=bool)
>>> mask[12, 12] = True
>>> aper = CircularAperture([(12.0, 12.0), (0.0, 12.0)], r=3.0)
>>> aperstats = ApertureStats(data, aper, mask=mask)
>>> for source_id, names in aperstats.decode_flags().items():
... print(source_id, names)
1 ['masked_pixels']
2 ['partial_overlap']
"""
decoded = decode_aperture_flags(self._array('flags'),
return_bit_values=return_bit_values)
return {int(id_): flags
for id_, flags in zip(self._array('id'), decoded,
strict=True)}
def _get_values(self, array):
"""
Get a 1D array of unmasked aperture-weighted values from the
input array.
An array with a single NaN is returned for completely-masked
sources.
"""
if self.isscalar:
array = (array,)
return [arr.compressed() if len(arr.compressed()) > 0
else np.array([np.nan]) for arr in array]
@cached_property
def _data_values_center(self):
"""
A 1D array of unmasked aperture-weighted data values using the
"center" method.
An array with a single NaN is returned for completely-masked
sources.
"""
return self._get_values(self.data_cutout)
[docs]
@cached_property
def moments(self):
"""
Spatial moments up to 3rd order of the source.
"""
gather = self._fast_gather
if gather is not None:
overlap = gather.overlap
zeros = np.zeros(self.n_positions)
mom = self._threaded_reduction(
batch_moments,
(gather.values, gather.local_x, gather.local_y),
gather.starts, gather.counts, per_source=(zeros, zeros))
# No-overlap sources have NaN moments (the mask-based path
# uses an all-NaN cutout); all-masked overlapping sources
# have zero moments (an all-zero cutout).
mom[~overlap] = np.nan
return mom
return np.array([_image_moments(arr, order=3)
for arr in self._moment_data_cutout])
[docs]
@cached_property
def moments_central(self):
"""
Central moments (translation invariant) of the source up to 3rd
order.
"""
cutout_centroid = self._array('cutout_centroid')
gather = self._fast_gather
if gather is not None:
cen_x = np.ascontiguousarray(cutout_centroid[:, 0])
cen_y = np.ascontiguousarray(cutout_centroid[:, 1])
mom = self._threaded_reduction(
batch_moments,
(gather.values, gather.local_x, gather.local_y),
gather.starts, gather.counts, per_source=(cen_x, cen_y))
# Empty sources (no overlap or fully masked) have a NaN
# centroid, so their central moments are NaN (matching the
# mask-based path).
mom[gather.counts == 0] = np.nan
return mom
return np.array([_image_moments(arr, center=(xcen_, ycen_), order=3)
for arr, xcen_, ycen_ in
zip(self._moment_data_cutout, cutout_centroid[:, 0],
cutout_centroid[:, 1], strict=True)])
[docs]
@cached_property
def cutout_centroid(self):
"""
The ``(x, y)`` coordinate, relative to the cutout data, of the
centroid within the aperture.
The centroid is computed as the center of mass of the unmasked
pixels within the aperture.
"""
moments = self._array('moments')
# Ignore divide-by-zero RuntimeWarning
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
y_centroid = moments[:, 1, 0] / moments[:, 0, 0]
x_centroid = moments[:, 0, 1] / moments[:, 0, 0]
return np.transpose((x_centroid, y_centroid))
[docs]
@cached_property
def centroid(self):
"""
The ``(x, y)`` coordinate of the centroid.
The centroid is computed as the center of mass of the unmasked
pixels within the aperture.
"""
origin = np.transpose((self.bbox_xmin, self.bbox_ymin))
return self.cutout_centroid + origin
[docs]
@cached_property
def x_centroid(self):
"""
The ``x`` coordinate of the centroid.
The centroid is computed as the center of mass of the unmasked
pixels within the aperture.
"""
return np.transpose(self._array('centroid'))[0]
[docs]
@cached_property
def y_centroid(self):
"""
The ``y`` coordinate of the centroid.
The centroid is computed as the center of mass of the unmasked
pixels within the aperture.
"""
return np.transpose(self._array('centroid'))[1]
[docs]
@cached_property
def sky_centroid(self):
"""
The sky coordinate of the centroid of the unmasked pixels within
the aperture, returned as a `~astropy.coordinates.SkyCoord`
object.
The output coordinate frame is the same as the input ``wcs``.
`None` if ``wcs`` is not input.
"""
if self._wcs is None:
return self._null_object
return self._wcs.pixel_to_world(self.x_centroid, self.y_centroid)
[docs]
@cached_property
def sky_centroid_icrs(self):
"""
The sky coordinate in the International Celestial Reference
System (ICRS) frame of the centroid of the unmasked pixels
within the aperture, returned as a
`~astropy.coordinates.SkyCoord` object.
`None` if ``wcs`` is not input.
"""
if self._wcs is None:
return self._null_object
return self.sky_centroid.icrs
@cached_property
def _bbox(self):
"""
The `~photutils.aperture.BoundingBox` of the aperture, always as
an iterable.
"""
apertures = self._pixel_aperture
if self.isscalar:
apertures = (apertures,)
return [aperture.bbox for aperture in apertures]
[docs]
@cached_property
def bbox(self):
"""
The `~photutils.aperture.BoundingBox` of the aperture.
Note that the aperture bounding box is calculated using the
exact size of the aperture, which may be slightly larger than
the aperture mask calculated using the "center" method.
"""
return self._bbox
@cached_property
def _bbox_bounds(self):
"""
The bounding box x and y minimum and maximum (inclusive)
bounds, as an ``(n_positions, 4)`` array.
When the aperture uses the default ``bbox`` implementation,
the bounds are taken from the aperture's vectorized integer
bounds, so no per-position `~photutils.aperture.BoundingBox`
objects are created. An aperture subclass that overrides
``bbox`` falls back to reading the per-position objects.
"""
aper = self._pixel_aperture
if type(aper).bbox is PixelAperture.bbox:
# Convert the exclusive upper bounds to inclusive
return aper._bbox_bounds - [0, 1, 0, 1]
bbox = self._array('bbox')
# The reshape preserves the (n_positions, 4) shape when there
# are zero positions
return np.array([(bbox_.ixmin, bbox_.ixmax - 1,
bbox_.iymin, bbox_.iymax - 1)
for bbox_ in bbox], dtype=int).reshape(-1, 4)
[docs]
@cached_property
def bbox_xmin(self):
"""
The minimum ``x``-pixel index of the bounding box.
"""
return np.transpose(self._bbox_bounds)[0]
[docs]
@cached_property
def bbox_xmax(self):
"""
The maximum ``x``-pixel index of the bounding box.
Note that this value is inclusive, unlike numpy slice indices.
"""
return np.transpose(self._bbox_bounds)[1]
[docs]
@cached_property
def bbox_ymin(self):
"""
The minimum ``y``-pixel index of the bounding box.
"""
return np.transpose(self._bbox_bounds)[2]
[docs]
@cached_property
def bbox_ymax(self):
"""
The maximum ``y``-pixel index of the bounding box.
Note that this value is inclusive, unlike numpy slice indices.
"""
return np.transpose(self._bbox_bounds)[3]
@cached_property
def _center_n_pixels(self):
"""
The number of unmasked pixels within each aperture using the
"center" mask method.
The result is a `~numpy.ndarray` of per-source pixel counts.
Sources with no unmasked pixels are set to NaN.
"""
gather = self._fast_gather
if gather is not None:
counts = gather.counts
n_pixels = counts.astype(float)
n_pixels[counts == 0] = np.nan
return n_pixels
n_pixels = np.array([np.sum(weight.filled(0.0))
for weight in self._weight_cutout_center])
n_pixels[self._all_masked] = np.nan
return n_pixels
@cached_property
def _sem(self):
"""
The standard error of the mean for each aperture.
The result is a `~numpy.ndarray` (without data units). It is
computed as ``s / sqrt(N)``, where ``s`` is the sample (``N -
1``) standard deviation and ``N`` is the number of unmasked
pixels within the aperture. Sources with fewer than two unmasked
pixels are set to NaN.
"""
stats = self._mean_var
if stats is not None:
var = stats[1]
else:
var = np.array([np.var(values)
for values in self._data_values_center])
n_pixels = self._center_n_pixels
sem = np.full(self.n_positions, np.nan)
mask = n_pixels >= 2
# var is the population (ddof=0) variance, so var / (N - 1)
# equals the squared standard error of the mean.
sem[mask] = np.sqrt(var[mask] / (n_pixels[mask] - 1.0))
return sem
[docs]
@cached_property
def center_aper_area(self):
"""
The total area of the unmasked pixels within the aperture using
the "center" aperture mask method.
"""
return self._center_n_pixels * (u.pix**2)
[docs]
@cached_property
def sum_aper_area(self):
"""
The total area of the unmasked pixels within the aperture using
the input ``sum_method`` aperture mask method.
"""
gather = self._fast_sum
if gather is not None:
area, overlap = gather.sum_area.copy(), gather.overlap
area[overlap & (area == 0)] = np.nan
return area << (u.pix**2)
areas = np.array([np.sum(weight.filled(0.0))
for weight in self._weight_cutout])
# NaN only when no sum-method pixel survives. The center-method
# ``_all_masked`` flag must not be used here: when ``sum_method``
# is not "center", unmasked boundary pixels can carry a nonzero
# fractional area even if every center-method pixel is masked
# (and their values then also contribute to ``sum``).
areas[areas == 0] = np.nan
return areas << (u.pix**2)
[docs]
@cached_property
def sum(self):
r"""
The sum of the unmasked ``data`` values within the aperture.
.. math::
F = \sum_{i \in A} I_i
where :math:`F` is ``sum``, :math:`I_i` is the
background-subtracted ``data``, and :math:`A` are the unmasked
pixels in the aperture.
Non-finite pixel values (NaN and inf) are excluded
(automatically masked).
"""
gather = self._fast_sum
if gather is not None:
result = gather.sum_aper.copy()
area = gather.sum_area
overlap = gather.overlap
# No sum-method survivors -> NaN (matches the mask-based
# path, which returns NaN for an all-masked aperture).
result[overlap & (area == 0)] = np.nan
if self._data_unit is not None:
result <<= self._data_unit
return result
if self.sum_method == 'center':
return self._finalize_value_stat(None, np.sum)
data_values = self._get_values(self.data_sum_cutout)
result = np.array([np.sum(arr) for arr in data_values])
if self._data_unit is not None:
result <<= self._data_unit
return result
[docs]
@cached_property
def sum_err(self):
r"""
The uncertainty of `sum`, propagated from the input ``error``
array.
Because `sum` is the sum of the pixel values weighted by their
aperture overlap fractions, ``sum_err`` is the quadrature
sum of the total errors over the unmasked pixels within the
aperture, with each pixel variance weighted by the squared
overlap fraction:
.. math::
\Delta F = \sqrt{\sum_{i \in A} w_i^2
\sigma_{\mathrm{tot}, i}^2}
where :math:`\Delta F` is the `sum_err`, :math:`w_i` are the
aperture overlap fractions (1 for pixels entirely within the
aperture), :math:`\sigma_{\mathrm{tot, i}}` are the pixel-wise
total errors (``error``), and :math:`A` are the unmasked pixels
in the aperture.
Pixel values that are masked in the input ``data``, including
any non-finite pixel values (NaN and inf) that are automatically
masked, are also masked in the error array.
"""
if self._error is None:
err = self._null_value
else:
gather = self._fast_sum
if gather is not None:
variance = gather.var_aper.copy()
area = gather.sum_area
overlap = gather.overlap
variance[overlap & (area == 0)] = np.nan
err = np.sqrt(variance)
else:
if self.sum_method == 'center':
variance = self._variance_cutout_center
else:
variance = self._variance_cutout
var_values = [arr.compressed() if len(arr.compressed()) > 0
else np.array([np.nan]) for arr in variance]
err = np.sqrt([np.sum(arr) for arr in var_values])
if self._data_unit is not None:
err <<= self._data_unit
return err
[docs]
@cached_property
def min(self):
"""
The minimum of the unmasked pixel values within the aperture.
"""
fast = None if self._order_stats is None else self._order_stats[0]
return self._finalize_value_stat(fast, np.min)
[docs]
@cached_property
def max(self):
"""
The maximum of the unmasked pixel values within the aperture.
"""
fast = None if self._order_stats is None else self._order_stats[1]
return self._finalize_value_stat(fast, np.max)
[docs]
@cached_property
def mean(self):
"""
The mean of the unmasked pixel values within the aperture.
"""
fast = None if self._mean_var is None else self._mean_var[0]
return self._finalize_value_stat(fast, np.mean)
[docs]
@cached_property
def mean_err(self):
r"""
The standard error of the `mean`.
``mean_err`` is the standard deviation of the sampling
distribution of the mean:
.. math::
\sigma_{\bar{x}} = \frac{s}{\sqrt{N}}
where :math:`s` is the sample standard deviation (computed with
``N - 1`` in the denominator) and :math:`N` is the number of
unmasked pixels within the aperture (`center_aper_area`).
Apertures with fewer than two unmasked pixels have an undefined
standard error and are set to NaN.
"""
result = self._sem.copy()
if self._data_unit is not None:
result <<= self._data_unit
return result
[docs]
@cached_property
def mode(self):
"""
The mode of the unmasked pixel values within the aperture.
The mode is estimated as ``(3 * median) - (2 * mean)``.
"""
return 3.0 * self.median - 2.0 * self.mean
@cached_property
def _variance(self):
"""
The variance of the unmasked pixel values within each aperture
as a plain `~numpy.ndarray` (without data units).
The population (``ddof=0``) variance is computed using the fast
Cython reduction when available and otherwise the mask-based
per-source code path. When ``self.ddof`` is nonzero, the result
is rescaled by ``N / (N - ddof)``, where ``N`` is the number of
unmasked pixels within the aperture. Apertures with ``N <=
ddof`` unmasked pixels are set to NaN.
"""
fast = None if self._mean_var is None else self._mean_var[1]
var = self._finalize_value_stat(fast, np.var, apply_unit=False)
if self.ddof == 0:
return var
n_pixels = self._center_n_pixels
result = np.full(self.n_positions, np.nan)
mask = n_pixels > self.ddof
result[mask] = (var[mask] * n_pixels[mask]
/ (n_pixels[mask] - self.ddof))
return result
[docs]
@cached_property
def std(self):
"""
The standard deviation of the unmasked pixel values within the
aperture.
The divisor used in the calculation is ``N - ddof``, where ``N``
is the number of unmasked pixels within the aperture and ``ddof``
is the value of the ``ddof`` keyword (default 0).
"""
result = np.sqrt(self._variance)
if self._data_unit is not None:
result <<= self._data_unit
return result
[docs]
@cached_property
def mad_std(self):
r"""
The standard deviation calculated using
the `median absolute deviation (MAD)
<https://en.wikipedia.org/wiki/Median_absolute_deviation>`_.
The standard deviation estimator is given by:
.. math::
\sigma \approx \frac{\textrm{MAD}}{\Phi^{-1}(3/4)}
\approx 1.4826 \ \textrm{MAD}
where :math:`\Phi^{-1}(P)` is the normal inverse cumulative
distribution function evaluated at probability :math:`P = 3/4`.
"""
fast = None if self._mad is None else self._mad * _MAD_STD_SCALE
return self._finalize_value_stat(fast, mad_std)
[docs]
@cached_property
def var(self):
"""
The variance of the unmasked pixel values within the aperture.
The divisor used in the calculation is ``N - ddof``, where ``N``
is the number of unmasked pixels within the aperture and ``ddof``
is the value of the ``ddof`` keyword (default 0).
"""
result = self._variance.copy()
if self._data_unit is not None:
result <<= self._data_unit**2
return result
[docs]
@cached_property
def biweight_location(self):
"""
The biweight location of the unmasked pixel values within the
aperture.
The tuning constant is fixed at ``c=6``, the default value
used by `astropy.stats.biweight_location`.
"""
fast = None if self._biweight is None else self._biweight[0]
return self._finalize_value_stat(fast, biweight_location)
[docs]
@cached_property
def biweight_midvariance(self):
"""
The biweight midvariance of the unmasked pixel values within the
aperture.
The tuning constant is fixed at ``c=9``, the default value
used by `astropy.stats.biweight_midvariance`.
"""
fast = None if self._biweight is None else self._biweight[1]
return self._finalize_value_stat(fast, biweight_midvariance,
square_unit=True)
[docs]
@cached_property
def inertia_tensor(self):
"""
The inertia tensor of the source for the rotation around its
center of mass.
"""
moments = self._array('moments_central')
mu_02 = moments[:, 0, 2]
mu_11 = -moments[:, 1, 1]
mu_20 = moments[:, 2, 0]
tensor = np.array([mu_02, mu_11, mu_11, mu_20]).swapaxes(0, 1)
return tensor.reshape((tensor.shape[0], 2, 2)) * u.pix**2
@cached_property
def _raw_covariance(self):
"""
The raw ``(N, 2, 2)`` covariance matrix of the 2D Gaussian
function that has the same normalized second-order moments as
the source, before any regularization.
This unregularized matrix is shared by `_covariance` (which
regularizes a copy) and `_singular_covariance_mask` (which tests
it for singularity). Callers that modify the matrix in place
must operate on a copy so the cached value is not corrupted.
"""
moments = self._array('moments_central')
# Ignore divide-by-zero RuntimeWarning
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
mu_norm = moments / moments[:, 0, 0][:, np.newaxis, np.newaxis]
covar = np.array([mu_norm[:, 0, 2], mu_norm[:, 1, 1],
mu_norm[:, 1, 1], mu_norm[:, 2, 0]]).swapaxes(0, 1)
return covar.reshape((covar.shape[0], 2, 2))
@cached_property
def _covariance(self):
"""
The covariance matrix of the 2D Gaussian function that has the
same second-order moments as the source, always as an iterable.
"""
# Copy so the regularization below does not mutate the cached
# raw covariance shared with `_singular_covariance_mask`.
covar = self._raw_covariance.copy()
# Regularize the covariance matrix for "infinitely" thin
# detections by incrementally increasing the diagonal elements
# by 1/12, the variance of a uniform distribution across a
# single pixel (the smallest second moment a resolved source can
# have given finite pixel size).
delta = 1.0 / 12
delta2 = delta**2
# Ignore RuntimeWarning from NaN values in covar
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
covar_det = np.linalg.det(covar)
covar_trace = covar[:, 0, 0] + covar[:, 1, 1]
# A valid covariance is positive semidefinite (det >= 0
# and trace >= 0). Any matrix that is not (e.g., from
# net-negative flux weighting) has an undefined shape and is
# set to NaN.
bad = (covar_det < 0) | (covar_trace < 0)
covar[bad] = np.nan
# Regularize "infinitely" thin detections by adding 1/12
# (delta) to each diagonal. A single bump is sufficient.
# For a positive semidefinite matrix the bumped determinant
# exceeds the raw determinant by delta times the trace plus
# delta squared. Since the raw determinant and trace are
# both non-negative, the result is at least delta squared,
# which equals the delta2 threshold, so it clears the
# threshold in a single step.
idx = np.where(covar_det < delta2)[0]
covar[idx, 0, 0] += delta
covar[idx, 1, 1] += delta
return covar
[docs]
@cached_property
def covariance(self):
"""
The covariance matrix of the 2D Gaussian function that has the
same second-order moments as the source.
"""
return self._covariance * (u.pix**2)
[docs]
@cached_property
def covariance_eigvals(self):
"""
The two eigenvalues of the `covariance` matrix in decreasing
order.
"""
eigvals = np.full((self.n_positions, 2), np.nan)
# np.linalg.eigvalsh requires that every element of a covariance
# matrix be finite, so select only the wholly finite matrices
idx = np.flatnonzero(np.isfinite(self._covariance).all(axis=(1, 2)))
eigvals[idx] = np.linalg.eigvalsh(self._covariance[idx])
# Check for negative variance
# (just in case covariance matrix is not positive semidefinite)
idx2 = np.unique(np.where(eigvals < 0)[0])
eigvals[idx2] = (np.nan, np.nan)
# Sort each eigenvalue pair in descending order
# (eigvalsh returns values in ascending order)
eigvals = np.fliplr(eigvals)
return eigvals * u.pix**2
[docs]
@cached_property
def semimajor_axis(self):
"""
The 1-sigma standard deviation along the semimajor axis of the
2D Gaussian function that has the same second-order central
moments as the source.
"""
eigvals = self._array('covariance_eigvals')
return np.sqrt(eigvals[:, 0])
[docs]
@cached_property
def semiminor_axis(self):
"""
The 1-sigma standard deviation along the semiminor axis of the
2D Gaussian function that has the same second-order central
moments as the source.
"""
eigvals = self._array('covariance_eigvals')
return np.sqrt(eigvals[:, 1])
[docs]
@cached_property
def fwhm(self):
r"""
The circularized full width at half maximum (FWHM) of the 2D
Gaussian function that has the same second-order central moments
as the source.
.. math::
\mathrm{FWHM} & = 2 \sqrt{2 \ln(2)} \sqrt{0.5 (a^2 + b^2)}
\\
& = 2 \sqrt{\ln(2) \ (a^2 + b^2)}
where :math:`a` and :math:`b` are the 1-sigma lengths of the
semimajor (`semimajor_axis`) and semiminor (`semiminor_axis`)
axes, respectively.
"""
return 2.0 * np.sqrt(np.log(2.0) * (self.semimajor_axis**2
+ self.semiminor_axis**2))
[docs]
@cached_property
def orientation(self):
"""
The angle between the ``x`` axis and the major axis of the 2D
Gaussian function that has the same second-order moments as the
source.
The angle increases in the counter-clockwise direction and is
in the range (-90, 90] degrees.
"""
covar = self._covariance
orient_radians = 0.5 * np.arctan2(2.0 * covar[:, 0, 1],
(covar[:, 0, 0] - covar[:, 1, 1]))
return np.rad2deg(orient_radians) * u.deg
[docs]
@cached_property
def eccentricity(self):
r"""
The eccentricity of the 2D Gaussian function that has the same
second-order moments as the source.
The eccentricity is the fraction of the distance along the
semimajor axis at which the focus lies.
.. math::
e = \sqrt{1 - \frac{b^2}{a^2}}
where :math:`a` and :math:`b` are the lengths of the semimajor
and semiminor axes, respectively.
"""
semimajor_var, semiminor_var = np.transpose(self.covariance_eigvals)
return np.sqrt(1.0 - (semiminor_var / semimajor_var))
[docs]
@cached_property
def elongation(self):
r"""
The ratio of the lengths of the semimajor and semiminor axes.
.. math::
\mathrm{elongation} = \frac{a}{b}
where :math:`a` and :math:`b` are the lengths of the semimajor
and semiminor axes, respectively.
"""
return self.semimajor_axis / self.semiminor_axis
[docs]
@cached_property
def ellipticity(self):
r"""
1.0 minus the ratio of the lengths of the semiminor and
semimajor axes (or 1.0 divided by the `elongation`, subtracted
from 1.0).
.. math::
\mathrm{ellipticity} = \frac{a - b}{a} = 1 - \frac{b}{a}
= 1 - \frac{1}{\mathrm{elongation}}
where :math:`a` and :math:`b` are the lengths of the semimajor
and semiminor axes, respectively.
"""
return 1.0 - (self.semiminor_axis / self.semimajor_axis)
[docs]
@cached_property
def covariance_xx(self):
r"""
The ``(0, 0)`` element of the `covariance` matrix, representing
:math:`\sigma_x^2`, in units of pixel**2.
"""
return self._covariance[:, 0, 0] * u.pix**2
[docs]
@cached_property
def covariance_yy(self):
r"""
The ``(1, 1)`` element of the `covariance` matrix, representing
:math:`\sigma_y^2`, in units of pixel**2.
"""
return self._covariance[:, 1, 1] * u.pix**2
[docs]
@cached_property
def covariance_xy(self):
r"""
The ``(0, 1)`` and ``(1, 0)`` elements of the `covariance`
matrix, representing :math:`\sigma_x \sigma_y`, in units of
pixel**2.
"""
return self._covariance[:, 0, 1] * u.pix**2
[docs]
@cached_property
def ellipse_cxx(self):
r"""
Coefficient for ``x**2`` in the generalized ellipse equation in
units of pixel**(-2).
The ellipse is defined as
.. math::
cxx (x - \bar{x})^2 + cxy (x - \bar{x}) (y - \bar{y}) +
cyy (y - \bar{y})^2 = R^2
where :math:`R` is a parameter which scales the ellipse (in
units of the axes lengths).
The isophotal limit of a source is well represented by :math:`R
\approx 3`.
"""
return ((np.cos(self.orientation) / self.semimajor_axis)**2
+ (np.sin(self.orientation) / self.semiminor_axis)**2)
[docs]
@cached_property
def ellipse_cyy(self):
r"""
Coefficient for ``y**2`` in the generalized ellipse equation in
units of pixel**(-2).
The ellipse is defined as
.. math::
cxx (x - \bar{x})^2 + cxy (x - \bar{x}) (y - \bar{y}) +
cyy (y - \bar{y})^2 = R^2
where :math:`R` is a parameter which scales the ellipse (in
units of the axes lengths).
The isophotal limit of a source is well represented by
:math:`R \approx 3`.
"""
return ((np.sin(self.orientation) / self.semimajor_axis)**2
+ (np.cos(self.orientation) / self.semiminor_axis)**2)
[docs]
@cached_property
def ellipse_cxy(self):
r"""
Coefficient for ``x * y`` in the generalized ellipse equation in
units of pixel**(-2).
The ellipse is defined as
.. math::
cxx (x - \bar{x})^2 + cxy (x - \bar{x}) (y - \bar{y}) +
cyy (y - \bar{y})^2 = R^2
where :math:`R` is a parameter which scales the ellipse (in
units of the axes lengths).
The isophotal limit of a source is well represented by :math:`R
\approx 3`.
"""
return (2.0 * np.cos(self.orientation) * np.sin(self.orientation)
* ((1.0 / self.semimajor_axis**2)
- (1.0 / self.semiminor_axis**2)))
[docs]
@cached_property
def gini(self):
r"""
The `Gini coefficient
<https://en.wikipedia.org/wiki/Gini_coefficient>`_ of the
unmasked pixel values within the aperture.
The Gini coefficient of the distribution of absolute flux values
is calculated using the prescription from `Lotz et al. 2004
<https://ui.adsabs.harvard.edu/abs/2004AJ....128..163L/abstract>`_
(Eq. 6) as:
.. math::
G = \frac{1}{\overline{|x|} \, n \, (n - 1)}
\sum^{n}_{i} (2i - n - 1) \left | x_i \right |
where :math:`\overline{|x|}` is the mean of the absolute value
of all pixel values :math:`x_i`. If the sum of all pixel values
is zero, the Gini coefficient is zero.
Negative pixel values are used via their absolute value. Invalid
values (NaN and inf) in the input are automatically excluded
from the calculation. If only a single finite pixel remains
after filtering, the Gini coefficient is 0.0.
"""
return self._finalize_value_stat(self._gini, gini_func,
apply_unit=False)
