stats.py

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
This module provides tools for calculating the properties of sources
defined by an Aperture.
"""

import functools
import inspect
import warnings
from copy import deepcopy

import astropy.units as u
import numpy as np
from astropy.nddata import NDData, StdDevUncertainty
from astropy.stats import (SigmaClip, biweight_location, biweight_midvariance,
                           mad_std)
from astropy.table import QTable
from astropy.utils import lazyproperty
from astropy.utils.exceptions import AstropyUserWarning

from photutils.aperture import Aperture, SkyAperture, region_to_aperture
from photutils.aperture.core import _aperture_metadata
from photutils.utils._misc import _get_meta
from photutils.utils._moments import _moments, _moments_central
from photutils.utils._quantity_helpers import process_quantities

__all__ = ['ApertureStats']


# default table columns for `to_table()` output
DEFAULT_COLUMNS = ['id', 'xcentroid', 'ycentroid', '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_sigma',
                   'semiminor_sigma', 'orientation', 'eccentricity']


def as_scalar(method):
    """
    Return a decorated method where it will always return a scalar value
    (instead of a length-1 tuple/list/array) if the class is scalar.

    Parameters
    ----------
    method : function
        The method to be decorated.

    Returns
    -------
    decorator : function
        The decorated method.
    """
    @functools.wraps(method)
    def _decorator(*args, **kwargs):
        result = method(*args, **kwargs)
        try:
            return (result[0] if args[0].isscalar and len(result) == 1
                    else result)
        except TypeError:  # if result has no len
            return result

    return _decorator


class ApertureStats:
    """
    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
    function.

    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 overlap of the aperture on
        the pixel grid. This method is used only for calculating the
        ``sum``, ``sum_error``, ``sum_aper_area``, ``data_sumcutout``,
        and ``error_sumcutout`` properties. All other properties use the
        "center" aperture mask method. Not all options are available for
        all aperture types. The following methods are available:

        * ``'exact'`` (default):
          The exact fractional overlap of the aperture and each pixel is
          calculated. The aperture weights will contain values between 0
          and 1.

        * ``'center'``:
          A pixel is considered to be entirely in or out of the aperture
          depending on whether its center is in or out of the aperture.
          The aperture weights will contain values only of 0 (out) and 1
          (in).

        * ``'subpixel'``:
          A pixel is divided into subpixels (see the ``subpixels``
          keyword), each of which are considered to be entirely in or
          out of the aperture depending on whether its center is in
          or out of the aperture. If ``subpixels=1``, this method is
          equivalent to ``'center'``. The aperture weights will contain
          values between 0 and 1.

    subpixels : int, optional
        For the ``'subpixel'`` method, resample pixels by this factor
        in each dimension. That is, each pixel is divided into
        ``subpixels**2`` subpixels. This keyword is ignored unless
        ``sum_method='subpixel'``.

    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.

    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`.

    Most source properties are calculated using the "center"
    aperture-mask method, which gives aperture weights of 0 or 1. This
    avoids the need to compute weighted statistics --- the ``data``
    pixel values are directly used.

    The input ``sum_method`` and ``subpixels`` keywords are used
    to determine the aperture-mask method when calculating the
    sum-related properties: ``sum``, ``sum_error``, ``sum_aper_area``,
    ``data_sumcutout``, and ``error_sumcutout``. The default is
    ``sum_method='exact'``, which produces exact aperture-weighted
    photometry.

    .. _SourceExtractor: https://sextractor.readthedocs.io/en/latest/

    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.xcentroid)  # doctest: +FLOAT_CMP
    149.99080259251238
    >>> print(aperstats.ycentroid)  # doctest: +FLOAT_CMP
    24.97484633000507
    >>> print(aperstats.centroid)  # doctest: +FLOAT_CMP
    [149.99080259  24.97484633]

    >>> print(aperstats.mean, aperstats.median)  # doctest: +FLOAT_CMP
    47.76300955780609 31.913789514433084

    >>> print(aperstats.std)  # doctest: +FLOAT_CMP
    39.193655383492974

    >>> print(aperstats.sum)  # doctest: +FLOAT_CMP
    9286.709206410273

    >>> print(aperstats.sum_aper_area) # doctest: +FLOAT_CMP
    201.0619298297468 pix2

    >>> # more than one aperture position
    >>> aper2 = CircularAperture(((150, 25), (90, 60)), 10)
    >>> aperstats2 = ApertureStats(data, aper2)
    >>> print(aperstats2.xcentroid)  # doctest: +FLOAT_CMP
    [149.98470724  89.97893946]
    >>> print(aperstats2.sum)  # doctest: +FLOAT_CMP
    [10177.62548482 36653.97704059]
    """

    def __init__(self, data, aperture, *, error=None, mask=None, wcs=None,
                 sigma_clip=None, sum_method='exact', subpixels=5,
                 local_bkg=None):

        if isinstance(data, NDData):
            data, error, mask, wcs = self._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 = self._validate_array(data, 'data', shape=False)
        self._data_unit = unit
        self._input_aperture = self._validate_aperture(aperture)
        aperture_meta = _aperture_metadata(aperture)  # use input aperture

        if isinstance(aperture, SkyAperture) and wcs is None:
            raise ValueError('A wcs is required when using a SkyAperture')

        # convert region to aperture if necessary
        if not isinstance(aperture, Aperture):
            aperture = region_to_aperture(aperture)
        self.aperture = aperture

        self._error = self._validate_array(error, 'error')
        self._mask = self._validate_array(mask, 'mask')
        self._wcs = wcs

        if sigma_clip is not None and not isinstance(sigma_clip, SigmaClip):
            raise TypeError('sigma_clip must be a SigmaClip instance')
        self.sigma_clip = sigma_clip

        self.sum_method = sum_method
        self.subpixels = subpixels

        self._local_bkg = np.zeros(self.n_apertures)  # no local bkg
        if local_bkg is not None:
            local_bkg = np.atleast_1d(local_bkg)
            if local_bkg.ndim != 1:
                raise ValueError('local_bkg must be a 1D array')

            n_local_bkg = len(local_bkg)
            if n_local_bkg not in (1, self.n_apertures):
                raise ValueError('local_bkg must be scalar or have the same '
                                 'length as the input aperture')
            local_bkg = np.broadcast_to(local_bkg, self.n_apertures)

            if np.any(~np.isfinite(local_bkg)):
                raise ValueError('local_bkg must not contain any non-finite '
                                 '(e.g., inf or NaN) values')
            self._local_bkg = local_bkg  # always an iterable

        self._ids = np.arange(self.n_apertures) + 1
        self.default_columns = DEFAULT_COLUMNS
        self.meta = _get_meta()
        self.meta.update(aperture_meta)

    @staticmethod
    def _unpack_nddata(data, error, mask, wcs):
        nddata_attr = {'error': error, 'mask': mask, 'wcs': wcs}
        for key, value in nddata_attr.items():
            if value is not None:
                warnings.warn(f'The {key!r} keyword will be ignored. Its '
                              'value is obtained from the input NDData '
                              'object.', AstropyUserWarning)

        mask = data.mask
        wcs = data.wcs

        if isinstance(data.uncertainty, StdDevUncertainty):
            if data.uncertainty.unit is None:
                error = data.uncertainty.array
            else:
                error = data.uncertainty.array * data.uncertainty.unit

        if data.unit is not None:
            data = u.Quantity(data.data, unit=data.unit)
        else:
            data = data.data

        return data, error, mask, wcs

    @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):
            raise TypeError('aperture must be an Aperture or Region object')
        return aperture

    def _validate_array(self, array, name, ndim=2, shape=True):
        if name == 'mask' and array is np.ma.nomask:
            array = None
        if array is not None:
            array = np.asanyarray(array)
            if array.ndim != ndim:
                raise ValueError(f'{name} must be a {ndim}D array.')
            if shape and array.shape != self._data.shape:
                raise ValueError(f'data and {name} must have the same shape.')
        return array

    @property
    def _lazyproperties(self):
        """
        A list of all class lazyproperties (even in superclasses).
        """
        def islazyproperty(obj):
            return isinstance(obj, lazyproperty)

        return [i[0] for i in inspect.getmembers(self.__class__,
                                                 predicate=islazyproperty)]

    @property
    def properties(self):
        """
        A sorted list of built-in source properties.
        """
        lazyproperties = [name for name in self._lazyproperties if not
                          name.startswith('_')]
        lazyproperties.sort()
        return lazyproperties

    def __getitem__(self, index):
        if self.isscalar:
            raise TypeError(f'A scalar {self.__class__.__name__!r} object '
                            'cannot be indexed')

        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',
                     'default_columns', 'meta')
        for attr in init_attr:
            setattr(newcls, attr, getattr(self, attr))

        # need to slice _aperture and _ids;
        # aperture determines isscalar (needed below)
        attrs = ('aperture', '_ids')
        for attr in attrs:
            setattr(newcls, attr, getattr(self, attr)[index])

        # slice evaluated lazyproperty objects
        keys = set(self.__dict__.keys()) & set(self._lazyproperties)
        keys.add('_local_bkg')  # iterable defined in __init__
        for key in keys:
            value = self.__dict__[key]

            # do not insert attributes that are always scalar (e.g.,
            # isscalar, n_apertures), 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__}>'
        with np.printoptions(threshold=25, edgeitems=5):
            fmt = [f'Length: {self.n_apertures}']
        return f'{cls_name}\n' + '\n'.join(fmt)

    def __repr__(self):
        return self.__str__()

    def __len__(self):
        if self.isscalar:
            raise TypeError(f'Scalar {self.__class__.__name__!r} object has '
                            'no len()')
        return self.n_apertures

    def __iter__(self):
        for item in range(len(self)):
            yield self.__getitem__(item)

    @lazyproperty
    def isscalar(self):
        """
        Whether the instance is scalar (e.g., a single aperture
        position).
        """
        return self._pixel_aperture.isscalar

    def copy(self):
        """
        Return a deep copy of this object.

        Returns
        -------
        result : `ApertureStats`
            A deep copy of this object.
        """
        return deepcopy(self)

    @lazyproperty
    def _null_object(self):
        """
        Return `None` values.
        """
        return np.array([None] * self.n_apertures)

    @lazyproperty
    def _null_value(self):
        """
        Return np.nan values.
        """
        values = np.empty(self.n_apertures)
        values.fill(np.nan)
        return values

    @property
    @as_scalar
    def id(self):
        """
        The aperture identification number(s).
        """
        return self._ids

    @property
    def ids(self):
        """
        The aperture identification number(s), always as an iterable
        `~numpy.ndarray`.
        """
        _ids = self._ids
        if self.isscalar:
            _ids = np.array((_ids,))
        return _ids

    def get_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.get_ids(id_num)

    def get_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.
        """
        for id_num in np.atleast_1d(id_nums):
            if id_num not in self.ids:
                raise ValueError(f'{id_num} is not a valid source ID number')

        sorter = np.argsort(self.id)
        indices = sorter[np.searchsorted(self.id, id_nums, sorter=sorter)]
        return self[indices]

    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.
        """
        if columns is None:
            table_columns = self.default_columns
        elif isinstance(columns, str):
            table_columns = [columns]
        else:
            table_columns = columns

        tbl = QTable()
        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,)

            tbl[column] = values
        return tbl

    @lazyproperty
    def n_apertures(self):
        """
        The number of positions in the input aperture.
        """
        if self.isscalar:
            return 1
        return len(self._pixel_aperture)

    @lazyproperty
    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

    @lazyproperty
    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

    @lazyproperty
    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

    @lazyproperty
    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

    @lazyproperty
    def _data_cutouts(self):
        """
        The local-background-subtracted unmasked data cutouts using the
        aperture bounding box, always as a 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):
        """
        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.

        Returns
        -------
        data, variance, mask, weights : list of `~numpy.ndarray`
            A list of cutout arrays for the data, variance, mask and weight
            arrays for each source (aperture position).
        """
        data_cutouts = []
        variance_cutouts = []
        mask_cutouts = []
        weight_cutouts = []
        overlaps = []

        for (data_cutout, apermask, slices) in zip(self._data_cutouts,
                                                   aperture_masks,
                                                   self._overlap_slices,
                                                   strict=True):

            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.
                data_mask = ~np.isfinite(data_cutout)
                if self._mask is not None:
                    data_mask |= self._mask[slc_large]

                overlap = True
                aperweight_cutout = apermask.data[slc_small]
                weight_cutout = aperweight_cutout * ~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 self.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 = self.sigma_clip(data_cutout_ma)

                    # define a mask of only the sigma-clipped pixels
                    sigclip_mask = data_sigclip.mask & ~mask_cutout
                    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 = self._error[slc_large]**2
                    variance_cutout = (variance * aperweight_cutout
                                       * ~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)

        # 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, strict=True))

    @lazyproperty
    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)

    @lazyproperty
    def _aperture_cutouts(self):
        """
        Aperture-weighted cutouts for the data, variance, total mask,
        and aperture weights using the input ``sum_method`` aperture
        mask method.
        """
        return self._make_aperture_cutouts(self._aperture_masks)

    @lazyproperty
    def _mask_cutout_center(self):
        """
        Boolean mask cutouts representing the total mask.

        The total mask is combination of the input ``mask``, non-finite
        ``data`` values, the cutout aperture mask using the "center"
        method, and the sigma-clip mask.
        """
        return list(zip(*self._aperture_cutouts_center, strict=True))[2]

    @lazyproperty
    def _mask_cutout(self):
        """
        Boolean mask cutouts representing the total mask.

        The total mask is 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]

    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_sumcutout`` mask.

        Units are not applied.
        """
        return [np.ma.masked_array(arr, mask=mask)
                for arr, mask in zip(array, self._mask_cutout, strict=True)]

    @lazyproperty
    @as_scalar
    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.
        """
        return self._make_masked_array_center(
            list(zip(*self._aperture_cutouts_center, strict=True))[0])

    @lazyproperty
    @as_scalar
    def data_sumcutout(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])

    @lazyproperty
    def _variance_cutout_center(self):
        """
        A 2D aperture-weighted variance cutout 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
        return self._make_masked_array_center(
            list(zip(*self._aperture_cutouts_center, strict=True))[1])

    @lazyproperty
    def _variance_cutout(self):
        """
        A 2D aperture-weighted variance cutout 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])

    @lazyproperty
    @as_scalar
    def error_sumcutout(self):
        """
        A 2D aperture-weighted error cutout 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 [np.sqrt(var) for var in self._variance_cutout]

    @lazyproperty
    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 of the aperture mask,
        pixels from the input ``mask``, non-finite ``data`` values (NaN
        and inf), and sigma-clipped pixels.
        """
        return self._make_masked_array_center(
            list(zip(*self._aperture_cutouts_center, strict=True))[3])

    @lazyproperty
    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 of 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])

    @lazyproperty
    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.data_cutout)  # self.data_cutout is a list
        if self.isscalar:
            data = (data,)

        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

    @lazyproperty
    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])

    @lazyproperty
    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 _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]

    @lazyproperty
    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)

    @lazyproperty
    @as_scalar
    def moments(self):
        """
        Spatial moments up to 3rd order of the source.
        """
        return np.array([_moments(arr, order=3) for arr in
                         self._moment_data_cutout])

    @lazyproperty
    @as_scalar
    def moments_central(self):
        """
        Central moments (translation invariant) of the source up to 3rd
        order.
        """
        cutout_centroid = self.cutout_centroid
        if self.isscalar:
            cutout_centroid = cutout_centroid[np.newaxis, :]
        return np.array([_moments_central(arr, center=(xcen_, ycen_), order=3)
                         for arr, xcen_, ycen_ in
                         zip(self._moment_data_cutout, cutout_centroid[:, 0],
                             cutout_centroid[:, 1], strict=True)])

    @lazyproperty
    @as_scalar
    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.moments
        if self.isscalar:
            moments = moments[np.newaxis, :]

        # ignore divide-by-zero RuntimeWarning
        with warnings.catch_warnings():
            warnings.simplefilter('ignore', RuntimeWarning)
            ycentroid = moments[:, 1, 0] / moments[:, 0, 0]
            xcentroid = moments[:, 0, 1] / moments[:, 0, 0]
        return np.transpose((xcentroid, ycentroid))

    @lazyproperty
    @as_scalar
    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

    @lazyproperty
    def _xcentroid(self):
        """
        The ``x`` coordinate of the centroid, always as an iterable.
        """
        xcentroid = np.transpose(self.centroid)[0]
        if self.isscalar:
            xcentroid = (xcentroid,)
        return xcentroid

    @lazyproperty
    @as_scalar
    def xcentroid(self):
        """
        The ``x`` coordinate of the centroid.

        The centroid is computed as the center of mass of the unmasked
        pixels within the aperture.
        """
        return self._xcentroid

    @lazyproperty
    def _ycentroid(self):
        """
        The ``y`` coordinate of the centroid, always as an iterable.
        """
        ycentroid = np.transpose(self.centroid)[1]
        if self.isscalar:
            ycentroid = (ycentroid,)
        return ycentroid

    @lazyproperty
    @as_scalar
    def ycentroid(self):
        """
        The ``y`` coordinate of the centroid.

        The centroid is computed as the center of mass of the unmasked
        pixels within the aperture.
        """
        return self._ycentroid

    @lazyproperty
    @as_scalar
    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.xcentroid, self.ycentroid)

    @lazyproperty
    @as_scalar
    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

    @lazyproperty
    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]

    @lazyproperty
    @as_scalar
    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" mode.
        """
        return self._bbox

    @lazyproperty
    @as_scalar
    def _bbox_bounds(self):
        """
        The bounding box x/y minimum and maximum bounds.
        """
        bbox = self.bbox
        if self.isscalar:
            bbox = (bbox,)
        return np.array([(bbox_.ixmin, bbox_.ixmax - 1,
                          bbox_.iymin, bbox_.iymax - 1)
                         for bbox_ in bbox])

    @lazyproperty
    @as_scalar
    def bbox_xmin(self):
        """
        The minimum ``x`` pixel index of the bounding box.
        """
        return np.transpose(self._bbox_bounds)[0]

    @lazyproperty
    @as_scalar
    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]

    @lazyproperty
    @as_scalar
    def bbox_ymin(self):
        """
        The minimum ``y`` pixel index of the bounding box.
        """
        return np.transpose(self._bbox_bounds)[2]

    @lazyproperty
    @as_scalar
    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]

    def _calculate_stats(self, stat_func, unit=None):
        """
        Apply the input ``stat_func`` to the 1D array of unmasked data
        values in the aperture.

        Units are applied if the input ``data`` has units.

        Parameters
        ----------
        stat_func : callable
            The callable to apply to the 1D `~numpy.ndarray` of unmasked
            data values.

        unit : `None` or `astropy.unit.Unit`, optional
            The unit to apply to the output data. This is used only
            if the input ``data`` has units. If `None` then the input
            ``data`` unit will be used.
        """
        result = np.array([stat_func(arr) for arr in self._data_values_center])
        if unit is None:
            unit = self._data_unit
        if unit is not None:
            result <<= unit
        return result

    @lazyproperty
    @as_scalar
    def center_aper_area(self):
        """
        The total area of the unmasked pixels within the aperture using
        the "center" aperture mask method.
        """
        areas = np.array([np.sum(weight.filled(0.0))
                          for weight in self._weight_cutout_center])
        areas[self._all_masked] = np.nan
        return areas << (u.pix**2)

    @lazyproperty
    @as_scalar
    def sum_aper_area(self):
        """
        The total area of the unmasked pixels within the aperture using
        the input ``sum_method`` aperture mask method.
        """
        areas = np.array([np.sum(weight.filled(0.0))
                          for weight in self._weight_cutout])
        areas[self._all_masked] = np.nan
        return areas << (u.pix**2)

    @lazyproperty
    @as_scalar
    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).
        """
        if self.sum_method == 'center':
            return self._calculate_stats(np.sum)

        data_values = self._get_values(self.data_sumcutout)
        result = np.array([np.sum(arr) for arr in data_values])
        if self._data_unit is not None:
            result <<= self._data_unit
        return result

    @lazyproperty
    @as_scalar
    def sum_err(self):
        r"""
        The uncertainty of `sum` , propagated from the input ``error``
        array.

        ``sum_err`` is the quadrature sum of the total errors over the
        unmasked pixels within the aperture:

        .. math::

            \Delta F = \sqrt{\sum_{i \in A} \sigma_{\mathrm{tot}, i}^2}

        where :math:`\Delta F` is the `sum_err`,
        :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:
            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

    @lazyproperty
    @as_scalar
    def min(self):
        """
        The minimum of the unmasked pixel values within the aperture.
        """
        return self._calculate_stats(np.min)

    @lazyproperty
    @as_scalar
    def max(self):
        """
        The maximum of the unmasked pixel values within the aperture.
        """
        return self._calculate_stats(np.max)

    @lazyproperty
    @as_scalar
    def mean(self):
        """
        The mean of the unmasked pixel values within the aperture.
        """
        return self._calculate_stats(np.mean)

    @lazyproperty
    @as_scalar
    def median(self):
        """
        The median of the unmasked pixel values within the aperture.
        """
        return self._calculate_stats(np.median)

    @lazyproperty
    @as_scalar
    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

    @lazyproperty
    @as_scalar
    def std(self):
        """
        The standard deviation of the unmasked pixel values within the
        aperture.
        """
        return self._calculate_stats(np.std)

    @lazyproperty
    @as_scalar
    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`.
        """
        return self._calculate_stats(mad_std)

    @lazyproperty
    @as_scalar
    def var(self):
        """
        The variance of the unmasked pixel values within the aperture.
        """
        unit = self._data_unit
        if unit is not None:
            unit **= 2
        return self._calculate_stats(np.var, unit=unit)

    @lazyproperty
    @as_scalar
    def biweight_location(self):
        """
        The biweight location of the unmasked pixel values within the
        aperture.

        See `astropy.stats.biweight_location`.
        """
        return self._calculate_stats(biweight_location)

    @lazyproperty
    @as_scalar
    def biweight_midvariance(self):
        """
        The biweight midvariance of the unmasked pixel values within the
        aperture.

        See `astropy.stats.biweight_midvariance`
        """
        unit = self._data_unit
        if unit is not None:
            unit **= 2
        return self._calculate_stats(biweight_midvariance, unit=unit)

    @lazyproperty
    @as_scalar
    def inertia_tensor(self):
        """
        The inertia tensor of the source for the rotation around its
        center of mass.
        """
        moments = self.moments_central
        if self.isscalar:
            moments = moments[np.newaxis, :]
        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

    @lazyproperty
    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.
        """
        moments = self.moments_central
        if self.isscalar:
            moments = moments[np.newaxis, :]
        # 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)
        covar = covar.reshape((covar.shape[0], 2, 2))

        # Modify the covariance matrix in the case of "infinitely" thin
        # detections. This follows SourceExtractor's prescription of
        # incrementally increasing the diagonal elements by 1/12.
        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)

            # covariance should be positive semidefinite
            idx = np.where(covar_det < 0)[0]
            covar[idx] = np.array([[np.nan, np.nan], [np.nan, np.nan]])

            idx = np.where(covar_det < delta2)[0]
            while idx.size > 0:
                covar[idx, 0, 0] += delta
                covar[idx, 1, 1] += delta
                covar_det = np.linalg.det(covar)
                idx = np.where(covar_det < delta2)[0]
        return covar

    @lazyproperty
    @as_scalar
    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)

    @lazyproperty
    @as_scalar
    def covariance_eigvals(self):
        """
        The two eigenvalues of the `covariance` matrix in decreasing
        order.
        """
        eigvals = np.empty((self.n_apertures, 2))
        eigvals.fill(np.nan)
        # np.linalg.eivals requires finite input values
        idx = np.unique(np.where(np.isfinite(self._covariance))[0])
        eigvals[idx] = np.linalg.eigvals(self._covariance[idx])

        # check for negative variance
        # (just in case covariance matrix is not positive semidefinite)
        idx2 = np.unique(np.where(eigvals < 0)[0])  # pragma: no cover
        eigvals[idx2] = (np.nan, np.nan)  # pragma: no cover

        # sort each eigenvalue pair in descending order
        eigvals.sort(axis=1)
        eigvals = np.fliplr(eigvals)

        return eigvals * u.pix**2

    @lazyproperty
    @as_scalar
    def semimajor_sigma(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.covariance_eigvals
        if self.isscalar:
            eigvals = eigvals[np.newaxis, :]
        # this matches SourceExtractor's A parameter
        return np.sqrt(eigvals[:, 0])

    @lazyproperty
    @as_scalar
    def semiminor_sigma(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.covariance_eigvals
        if self.isscalar:
            eigvals = eigvals[np.newaxis, :]
        # this matches SourceExtractor's B parameter
        return np.sqrt(eigvals[:, 1])

    @lazyproperty
    @as_scalar
    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_sigma`) and semiminor (`semiminor_sigma`)
        axes, respectively.
        """
        return 2.0 * np.sqrt(np.log(2.0) * (self.semimajor_sigma**2
                                            + self.semiminor_sigma**2))

    @lazyproperty
    @as_scalar
    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.
        """
        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

    @lazyproperty
    @as_scalar
    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))

    @lazyproperty
    @as_scalar
    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_sigma / self.semiminor_sigma

    @lazyproperty
    @as_scalar
    def ellipticity(self):
        r"""
        1.0 minus the ratio of the lengths of the semimajor and
        semiminor axes (or 1.0 minus the `elongation`).

        .. math::

            \mathrm{ellipticity} = 1 - \frac{b}{a}

        where :math:`a` and :math:`b` are the lengths of the semimajor
        and semiminor axes, respectively.
        """
        return 1.0 - (self.semiminor_sigma / self.semimajor_sigma)

    @lazyproperty
    @as_scalar
    def covar_sigx2(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

    @lazyproperty
    @as_scalar
    def covar_sigy2(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

    @lazyproperty
    @as_scalar
    def covar_sigxy(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

    @lazyproperty
    @as_scalar
    def 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).

        `SourceExtractor`_ reports that the isophotal limit of a source
        is well represented by :math:`R \approx 3`.
        """
        return ((np.cos(self.orientation) / self.semimajor_sigma)**2
                + (np.sin(self.orientation) / self.semiminor_sigma)**2)

    @lazyproperty
    @as_scalar
    def 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).

        `SourceExtractor`_ reports that the isophotal limit of a source
        is well represented by :math:`R \approx 3`.
        """
        return ((np.sin(self.orientation) / self.semimajor_sigma)**2
                + (np.cos(self.orientation) / self.semiminor_sigma)**2)

    @lazyproperty
    @as_scalar
    def 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).

        `SourceExtractor`_ reports that 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_sigma**2)
                   - (1.0 / self.semiminor_sigma**2)))

    @lazyproperty
    @as_scalar
    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 is calculated using the prescription from
        `Lotz et al. 2004
        <https://ui.adsabs.harvard.edu/abs/2004AJ....128..163L/abstract>`_
        as:

        .. math::

            G = \frac{1}{\left | \bar{x} \right | n (n - 1)}
                \sum^{n}_{i} (2i - n - 1) \left | x_i \right |

        where :math:`\bar{x}` is the mean over pixel values :math:`x_i`
        within the aperture.

        The Gini coefficient is a way of measuring the inequality in a
        given set of values. In the context of galaxy morphology, it
        measures how the light of a galaxy image is distributed among
        its pixels. A Gini coefficient value of 0 corresponds to a
        galaxy image with the light evenly distributed over all pixels
        while a Gini coefficient value of 1 represents a galaxy image
        with all its light concentrated in just one pixel.
        """
        gini = []
        for arr in self._data_values_center:
            if np.all(np.isnan(arr)):
                gini.append(np.nan)
                continue
            npix = np.size(arr)
            normalization = np.abs(np.mean(arr)) * npix * (npix - 1)
            kernel = ((2.0 * np.arange(1, npix + 1) - npix - 1)
                      * np.abs(np.sort(arr)))
            gini.append(np.sum(kernel) / normalization)
        return np.array(gini)