# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Tools for centroiding sources.
"""
import inspect
import warnings
import numpy as np
from astropy.nddata import overlap_slices
from astropy.utils.exceptions import AstropyUserWarning
from photutils.centroids._utils import _process_data_mask
from photutils.utils._deprecation import (deprecated_positional_kwargs,
deprecated_renamed_argument)
from photutils.utils._parameters import as_pair
from photutils.utils._quantity_helpers import process_quantities
from photutils.utils._repr import make_repr
from photutils.utils._round import round_half_away
__all__ = ['CentroidQuadratic', 'centroid_com', 'centroid_quadratic',
'centroid_sources']
[docs]
@deprecated_positional_kwargs(since='3.0', until='4.0')
def centroid_com(data, mask=None):
"""
Calculate the centroid of an array as the flux-weighted
center of mass derived from `image moments
<https://en.wikipedia.org/wiki/Image_moment>`_.
Non-finite values (e.g., NaN or inf) in the ``data`` array are
automatically masked. The final mask is a logical OR combination
of the input ``mask``, the automatically generated mask for
non-finite values, and the mask of the input ``data`` if it is a
`~numpy.ma.MaskedArray`. The centroid is calculated using only the
unmasked data values.
Parameters
----------
data : array_like
The input n-dimensional array. ``data`` can be a
`~numpy.ma.MaskedArray`. The image should be a
background-subtracted cutout image containing a single
source. The source should be significantly stronger than the
background noise. If the data contains nearly equal positive and
negative values (i.e., the sum is close to zero), the centroid
calculation will be numerically unstable and may produce
undefined results that fall outside the array bounds.
mask : bool `~numpy.ndarray`, optional
A boolean mask, with the same shape as ``data``, where a `True`
value indicates the corresponding element of ``data`` is masked.
If ``data`` is a `~numpy.ma.MaskedArray`, its mask will be
combined (using bitwise OR) with the input ``mask``.
Returns
-------
centroid : `~numpy.ndarray`
The coordinates of the centroid in pixel order (e.g., ``(x, y)``
or ``(x, y, z)``), not numpy axis order. If the absolute value
of the sum of the (unmasked) data is smaller than 1e-30 (i.e.,
consistent with zero), then a `~numpy.ndarray` of NaN values
will be returned. If the sum is close to zero, the centroid may
be poorly defined and fall outside the array bounds.
Notes
-----
The centroid is calculated as:
.. math::
x_c = \\frac{\\sum x_i I_i}{\\sum I_i}, \\quad
y_c = \\frac{\\sum y_i I_i}{\\sum I_i}
where :math:`I_i` is the intensity at pixel :math:`(x_i, y_i)`.
Examples
--------
>>> import numpy as np
>>> from photutils.datasets import make_4gaussians_image
>>> from photutils.centroids import centroid_com
>>> data = make_4gaussians_image()
>>> data -= np.median(data[0:30, 0:125])
>>> data = data[40:80, 70:110]
>>> x1, y1 = centroid_com(data)
>>> print(np.array((x1, y1)))
[19.9796724 20.00992593]
.. plot::
import matplotlib.pyplot as plt
import numpy as np
from photutils.centroids import centroid_com
from photutils.datasets import make_4gaussians_image
data = make_4gaussians_image()
data -= np.median(data[0:30, 0:125])
data = data[40:80, 70:110]
xycen = centroid_com(data)
fig, ax = plt.subplots(figsize=(8, 8))
ax.imshow(data, origin='lower')
ax.scatter(*xycen, color='red', marker='+', s=100, label='Centroid')
ax.legend()
"""
(data,), _ = process_quantities((data,), ('data',))
data = _process_data_mask(data, mask, ndim=None, fill_value=0.0)
total = np.sum(data)
if abs(total) < 1.e-30:
return np.full(data.ndim, np.nan)
indices = np.ogrid[tuple(slice(0, i) for i in data.shape)]
# Output array is reversed to give (x, y) order (e.g., for 2D data)
return np.array([np.sum(indices[axis] * data) / total
for axis in range(data.ndim)])[::-1]
[docs]
@deprecated_positional_kwargs(since='3.0', until='4.0')
@deprecated_renamed_argument('xpeak', None, '3.0', until='4.0')
@deprecated_renamed_argument('ypeak', None, '3.0', until='4.0')
@deprecated_renamed_argument('search_boxsize', None, '3.0', until='4.0')
def centroid_quadratic(data, mask=None, fit_boxsize=5, xpeak=None,
ypeak=None, search_boxsize=None):
"""
Calculate the centroid of a 2D array by fitting a 2D quadratic
polynomial.
Non-finite values (e.g., NaN or inf) in the ``data`` array are
automatically masked. The final mask is a logical OR combination
of the input ``mask``, the automatically generated mask for
non-finite values, and the mask of the input ``data`` if it is a
`~numpy.ma.MaskedArray`. The centroid is calculated using only the
unmasked data values.
A second degree 2D polynomial is fit within a small region of the
data defined by ``fit_boxsize`` to calculate the centroid position.
The initial center of the fitting box can be specified using the
``xpeak`` and ``ypeak`` keywords. If both ``xpeak`` and ``ypeak``
are `None`, then the box will be centered at the position of the
maximum value in the input ``data``.
If ``xpeak`` and ``ypeak`` are specified, the ``search_boxsize``
optional keyword can be used to further refine the initial center of
the fitting box by searching for the position of the maximum pixel
within a box of size ``search_boxsize``.
`Vakili & Hogg (2016) <https://arxiv.org/abs/1610.05873>`_
demonstrate that 2D quadratic centroiding comes very
close to saturating the `Cramér-Rao lower bound
<https://en.wikipedia.org/wiki/Cram%C3%A9r%E2%80%93Rao_bound>`_ in a
wide range of conditions.
Parameters
----------
data : 2D array_like
The 2D image data. ``data`` can be a `~numpy.ma.MaskedArray`.
The image should be a background-subtracted cutout image
containing a single source.
mask : 2D bool `~numpy.ndarray`, 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 calculations. If ``data`` is
a `~numpy.ma.MaskedArray`, its mask will be combined (using
bitwise OR) with the input ``mask``.
fit_boxsize : int or tuple of int, optional
The size (in pixels) of the box used to define the fitting
region. If ``fit_boxsize`` has two elements, they must be in
``(ny, nx)`` order. If ``fit_boxsize`` is a scalar then a square
box of size ``fit_boxsize`` will be used. ``fit_boxsize`` must
have odd values for both axes.
xpeak, ypeak : float or `None`, optional
The initial guess of the position of the centroid. If either
``xpeak`` or ``ypeak`` is `None` then the position of the
maximum value in the input ``data`` will be used as the initial
guess.
.. deprecated:: 3.0
The ``xpeak`` and ``ypeak`` keywords are deprecated
and will be removed in a future version. Use
`~photutils.centroids.centroid_sources` to centroid sources
at specific positions.
search_boxsize : int or tuple of int, optional
The size (in pixels) of the box used to search for the maximum
pixel value if ``xpeak`` and ``ypeak`` are both specified. If
``search_boxsize`` has two elements, they must be in ``(ny,
nx)`` order. If ``search_boxsize`` is a scalar then a square
box of size ``search_boxsize`` will be used. ``search_boxsize``
must have odd values for both axes. This parameter is ignored
if either ``xpeak`` or ``ypeak`` is `None`. In that case, the
entire array is searched for the maximum value.
.. deprecated:: 3.0
The ``search_boxsize`` keyword is deprecated
and will be removed in a future version. Use
`~photutils.centroids.centroid_sources` to centroid sources
at specific positions.
Returns
-------
centroid : `~numpy.ndarray`
The ``x, y`` coordinates of the centroid.
Notes
-----
Use ``fit_boxsize = (3, 3)`` to match the work of `Vakili &
Hogg (2016) <https://arxiv.org/abs/1610.05873>`_ for their 2D
second-order polynomial centroiding method.
Because this centroid is based on fitting data, it can fail for many
reasons, returning (np.nan, np.nan):
* quadratic fit failed
* quadratic fit does not have a maximum
* quadratic fit maximum falls outside image
* not enough unmasked data points (6 are required)
A `ValueError` is raised if all data values are masked or
non-finite.
Also note that a fit is not performed if the maximum data value is
at the edge of the data. In this case, the position of the maximum
pixel will be returned.
References
----------
.. [1] Vakili and Hogg 2016, "Do fast stellar centroiding methods
saturate the Cramér-Rao lower bound?", `arXiv:1610.05873
<https://arxiv.org/abs/1610.05873>`_
Examples
--------
>>> import numpy as np
>>> from photutils.datasets import make_4gaussians_image
>>> from photutils.centroids import centroid_quadratic
>>> data = make_4gaussians_image()
>>> data -= np.median(data[0:30, 0:125])
>>> data = data[40:80, 70:110]
>>> x1, y1 = centroid_quadratic(data)
>>> print(np.array((x1, y1)))
[19.94009505 20.06884997]
.. plot::
import matplotlib.pyplot as plt
import numpy as np
from photutils.centroids import centroid_quadratic
from photutils.datasets import make_4gaussians_image
data = make_4gaussians_image()
data -= np.median(data[0:30, 0:125])
data = data[40:80, 70:110]
xycen = centroid_quadratic(data)
fig, ax = plt.subplots(figsize=(8, 8))
ax.imshow(data, origin='lower')
ax.scatter(*xycen, color='red', marker='+', s=100, label='Centroid')
ax.legend()
"""
(data,), _ = process_quantities((data,), ('data',))
if ((xpeak is None and ypeak is not None)
or (xpeak is not None and ypeak is None)):
msg = 'xpeak and ypeak must both be input or "None"'
raise ValueError(msg)
data = _process_data_mask(data, mask)
ny, nx = data.shape
if not np.any(np.isfinite(data)):
msg = 'All data values are masked or non-finite'
raise ValueError(msg)
fit_boxsize = as_pair('fit_boxsize', fit_boxsize, lower_bound=(0, 0),
upper_bound=data.shape, check_odd=True)
if np.prod(fit_boxsize) < 6:
msg = ('fit_boxsize is too small. 6 values are required to fit a '
'2D quadratic polynomial.')
raise ValueError(msg)
if xpeak is not None and ((xpeak < 0) or (xpeak > data.shape[1] - 1)):
msg = 'xpeak is outside the input data'
raise ValueError(msg)
if ypeak is not None and ((ypeak < 0) or (ypeak > data.shape[0] - 1)):
msg = 'ypeak is outside the input data'
raise ValueError(msg)
if xpeak is None or ypeak is None:
yidx, xidx = np.unravel_index(np.nanargmax(data), data.shape)
else:
xidx = round_half_away(xpeak)
yidx = round_half_away(ypeak)
if search_boxsize is not None:
search_boxsize = as_pair('search_boxsize', search_boxsize,
lower_bound=(0, 0),
upper_bound=data.shape, check_odd=True)
slc_data, _ = overlap_slices(data.shape, search_boxsize,
(yidx, xidx), mode='trim')
cutout = data[slc_data]
yidx, xidx = np.unravel_index(np.nanargmax(cutout), cutout.shape)
xidx += slc_data[1].start
yidx += slc_data[0].start
# Return the position of the maximum if it is at the edge of the
# data
if xidx in (0, nx - 1) or yidx in (0, ny - 1):
msg = ('maximum value is at the edge of the data and its '
'position was returned; no quadratic fit was performed')
warnings.warn(msg, AstropyUserWarning)
return np.array((xidx, yidx), dtype=float)
# Extract the fitting region
slc_data, _ = overlap_slices(data.shape, fit_boxsize, (yidx, xidx),
mode='trim')
xidx0, xidx1 = (slc_data[1].start, slc_data[1].stop)
yidx0, yidx1 = (slc_data[0].start, slc_data[0].stop)
# Shift the fitting box if it was clipped by the data edge
if (xidx1 - xidx0) < fit_boxsize[1]:
if xidx0 == 0:
xidx1 = min(nx, xidx0 + fit_boxsize[1])
if xidx1 == nx:
xidx0 = max(0, xidx1 - fit_boxsize[1])
if (yidx1 - yidx0) < fit_boxsize[0]:
if yidx0 == 0:
yidx1 = min(ny, yidx0 + fit_boxsize[0])
if yidx1 == ny:
yidx0 = max(0, yidx1 - fit_boxsize[0])
cutout = data[yidx0:yidx1, xidx0:xidx1].ravel()
if np.count_nonzero(~np.isnan(cutout)) < 6:
msg = ('at least 6 unmasked data points are required to '
'perform a 2D quadratic fit')
warnings.warn(msg, AstropyUserWarning)
return np.array((np.nan, np.nan))
# Fit a 2D quadratic polynomial to the fitting region. The fit
# coordinates are centered on the peak pixel to keep the design
# matrix well conditioned. With absolute coordinates the condition
# number grows as ~coordinate**4 and the fit fails for sources at
# large pixel coordinates (e.g., in large mosaic images).
xi = np.arange(xidx0, xidx1) - xidx
yi = np.arange(yidx0, yidx1) - yidx
x, y = np.meshgrid(xi, yi)
x = x.ravel()
y = y.ravel()
# Pre-allocate coefficient matrix for optimization
coeff_matrix = np.empty((x.size, 6), dtype=float)
coeff_matrix[:, 0] = 1
coeff_matrix[:, 1] = x
coeff_matrix[:, 2] = y
coeff_matrix[:, 3] = x * y
coeff_matrix[:, 4] = x * x
coeff_matrix[:, 5] = y * y
# Include only finite values in the fit.
finite_mask = np.isfinite(cutout)
if not np.all(finite_mask):
coeff_matrix = coeff_matrix[finite_mask]
cutout = cutout[finite_mask]
try:
c = np.linalg.lstsq(coeff_matrix, cutout, rcond=None)[0]
except np.linalg.LinAlgError:
msg = 'quadratic fit failed'
warnings.warn(msg, AstropyUserWarning)
return np.array((np.nan, np.nan))
# Analytically find the maximum of the polynomial
_, c10, c01, c11, c20, c02 = c
det = 4 * c20 * c02 - c11**2
# If the determinant is <= 0, the surface has a saddle point. If
# the determinant is > 0, the surface has a minimum or maximum. The
# curvature is negative (maximum) if c20 < 0 and c02 < 0. However,
# if det > 0, then 4 * c20 * c02 > c11**2 >= 0, so c20 and c02 must
# have the same sign. Therefore, we only need to check if c20 > 0
# (or c02 > 0) to determine if the surface has a minimum.
if det <= 0 or c20 > 0:
msg = 'quadratic fit does not have a maximum'
warnings.warn(msg, AstropyUserWarning)
return np.array((np.nan, np.nan))
# Add back the peak-pixel offset to convert the analytic maximum
# from fit coordinates to data coordinates
xm = (c01 * c11 - 2.0 * c02 * c10) / det + xidx
ym = (c10 * c11 - 2.0 * c20 * c01) / det + yidx
if 0.0 < xm < (nx - 1.0) and 0.0 < ym < (ny - 1.0):
xycen = np.array((xm, ym), dtype=float)
else:
msg = 'quadratic polynomial maximum value falls outside of the image'
warnings.warn(msg, AstropyUserWarning)
return np.array((np.nan, np.nan))
return xycen
[docs]
class CentroidQuadratic:
"""
Class to calculate the centroid of a 2D array by fitting a 2D
quadratic polynomial.
This class provides a callable interface to the
`~photutils.centroids.centroid_quadratic` function, allowing a
centroid function with specific fit parameters to be defined and
reused. This is useful, for example, when using a customized
centroid function with `~photutils.centroids.centroid_sources`.
Parameters
----------
fit_boxsize : int or tuple of int, optional
The size (in pixels) of the box used to define the fitting
region. If ``fit_boxsize`` has two elements, they must be in
``(ny, nx)`` order. If ``fit_boxsize`` is a scalar then a square
box of size ``fit_boxsize`` will be used. ``fit_boxsize`` must
have odd values for both axes.
Examples
--------
>>> import numpy as np
>>> from photutils.datasets import make_4gaussians_image
>>> from photutils.centroids import CentroidQuadratic
>>> data = make_4gaussians_image()
>>> data -= np.median(data[0:30, 0:125])
>>> data = data[40:80, 70:110]
>>> centroid_func = CentroidQuadratic(fit_boxsize=5)
>>> x1, y1 = centroid_func(data)
>>> print(np.array((x1, y1)))
[19.94009505 20.06884997]
Using with `~photutils.centroids.centroid_sources`::
>>> from photutils.centroids import centroid_sources
>>> data = make_4gaussians_image()
>>> data -= np.median(data[0:30, 0:125])
>>> x_init = (25, 91, 151, 160)
>>> y_init = (40, 61, 24, 71)
>>> centroid_func = CentroidQuadratic(fit_boxsize=3)
>>> x, y = centroid_sources(data, x_init, y_init, box_size=25,
... centroid_func=centroid_func)
"""
def __init__(self, *, fit_boxsize=5):
self.fit_boxsize = fit_boxsize
def __repr__(self):
return make_repr(self, ['fit_boxsize'])
def __str__(self):
return make_repr(self, ['fit_boxsize'], long=True)
[docs]
def __call__(self, data, *, mask=None):
"""
Calculate the centroid.
Non-finite values (e.g., NaN or inf) in the ``data`` array
are automatically masked. The automatically masked values are
combined (using bitwise OR) with the input ``mask``. If ``data``
is a `~numpy.ma.MaskedArray`, its mask will also be combined
(using bitwise OR) with the input ``mask``.
Parameters
----------
data : 2D array_like
The 2D image data. ``data`` can be a
`~numpy.ma.MaskedArray`. The image should be a
background-subtracted cutout image containing a single
source.
mask : 2D bool `~numpy.ndarray`, optional
A boolean mask, with the same shape as ``data``, where a
`True` value indicates the corresponding element of ``data``
is masked. If ``data`` is a `~numpy.ma.MaskedArray`, its
mask will be combined (using bitwise OR) with the input
``mask``. Masked data are excluded from calculations.
Returns
-------
centroid : `~numpy.ndarray`
The ``x, y`` coordinates of the centroid.
Notes
-----
Unlike `~photutils.centroids.centroid_1dg` and
`~photutils.centroids.centroid_2dg`, this method does not
support an error array.
"""
kwargs = {'mask': mask,
'fit_boxsize': self.fit_boxsize,
}
return centroid_quadratic(data, **kwargs)
[docs]
@deprecated_positional_kwargs(since='3.0', until='4.0')
def centroid_sources(data, xpos, ypos, box_size=11, footprint=None,
mask=None, centroid_func=centroid_com, **kwargs):
"""
Calculate the centroid of sources at the defined positions in a 2D
array using a specified centroid function.
A cutout image centered on each input position will be used to
calculate the centroid position. The cutout image is defined either
using the ``box_size`` or ``footprint`` keyword. The ``footprint``
keyword can be used to create a non-rectangular cutout image.
Masks and non-finite values are handled by the input
``centroid_func``. When using a centroid function provided by
Photutils, non-finite values (e.g., NaN or inf) in the ``data``
array are automatically masked. The ``centroid_1dg`` and
``centroid_2dg`` functions also automatically mask any pixels with
non-finite ``error`` array values. The final mask is a logical OR
combination of the input ``mask``, the automatically generated
mask(s) for non-finite values, and the mask of the input ``data`` if
it is a `~numpy.ma.MaskedArray`. The centroid is calculated using
only the unmasked data values.
Parameters
----------
data : 2D array_like
The 2D image data. ``data`` can be a `~numpy.ma.MaskedArray`.
The image should be background-subtracted.
xpos, ypos : float or array_like of float
The initial ``x`` and ``y`` pixel position(s) of the center
position. A cutout image centered on this position will be used
to calculate the centroid.
box_size : int or array_like of int, optional
The size of the cutout image along each axis. If ``box_size`` is
a number, then a square cutout of ``box_size`` will be created.
If ``box_size`` has two elements, they must be in ``(ny, nx)``
order. ``box_size`` must have odd values for both axes. Either
``box_size`` or ``footprint`` must be defined. If they are both
defined, then ``footprint`` overrides ``box_size``.
footprint : bool `~numpy.ndarray`, optional
A 2D boolean array where `True` values describe the local
footprint region to cutout. ``footprint`` can be used to create
a non-rectangular cutout image, in which case the input ``xpos``
and ``ypos`` represent the center of the minimal bounding box
for the input ``footprint``. ``box_size=(n, m)`` is equivalent
to ``footprint=np.ones((n, m))``. Either ``box_size`` or
``footprint`` must be defined. If they are both defined, then
``footprint`` overrides ``box_size``. The same ``footprint`` is
used for all sources.
mask : 2D bool `~numpy.ndarray`, optional
A 2D boolean array with the same shape as ``data``, where a
`True` value indicates the corresponding element of ``data`` is
masked. If ``data`` is a `~numpy.ma.MaskedArray`, its mask will
be combined (using bitwise OR) with the input ``mask``.
centroid_func : callable, optional
A callable object (e.g., function or class) that is used to
calculate the centroid of a 2D array. The ``centroid_func``
must accept a 2D `~numpy.ndarray`, have a ``mask`` keyword and
optionally an ``error`` keyword. A callable whose signature
accepts arbitrary keyword arguments (``**kwargs``) is assumed to
handle a ``mask`` keyword. The callable object must return two
scalar values representing the (x, y) centroid. The default is
`~photutils.centroids.centroid_com`.
**kwargs : dict, optional
Any additional keyword arguments accepted by the
``centroid_func``. A `TypeError` is raised for keyword arguments
not accepted by the ``centroid_func``.
Returns
-------
xcentroid, ycentroid : `~numpy.ndarray`
The ``x`` and ``y`` pixel position(s) of the centroids. NaNs
will be returned where the centroid failed. This is usually due
to a ``box_size`` that is too small when using a fitting-based
centroid function (e.g., `centroid_1dg`, `centroid_2dg`, or
`centroid_quadratic`).
Examples
--------
>>> import numpy as np
>>> from photutils.centroids import centroid_2dg, centroid_sources
>>> from photutils.datasets import make_4gaussians_image
>>> data = make_4gaussians_image()
>>> data -= np.median(data[0:30, 0:125])
>>> x_init = (25, 91, 151, 160)
>>> y_init = (40, 61, 24, 71)
>>> x, y = centroid_sources(data, x_init, y_init, box_size=25,
... centroid_func=centroid_2dg)
>>> print(x)
[ 24.96807828 89.98684636 149.96545721 160.18810915]
>>> print(y)
[40.03657613 60.01836631 24.96777946 69.80208702]
.. plot::
import matplotlib.pyplot as plt
import numpy as np
from photutils.centroids import centroid_2dg, centroid_sources
from photutils.datasets import make_4gaussians_image
data = make_4gaussians_image()
data -= np.median(data[0:30, 0:125])
x_init = (25, 91, 151, 160)
y_init = (40, 61, 24, 71)
x, y = centroid_sources(data, x_init, y_init, box_size=25,
centroid_func=centroid_2dg)
fig, ax = plt.subplots(figsize=(8, 4))
ax.imshow(data, origin='lower')
ax.scatter(x, y, marker='+', s=80, color='red', label='Centroids')
ax.legend()
fig.tight_layout()
"""
if np.ndim(data) != 2:
msg = 'data must be a 2D array'
raise ValueError(msg)
xpos = np.atleast_1d(xpos)
ypos = np.atleast_1d(ypos)
if xpos.ndim != 1:
msg = 'xpos must be a 1D array'
raise ValueError(msg)
if ypos.ndim != 1:
msg = 'ypos must be a 1D array'
raise ValueError(msg)
if len(xpos) != len(ypos):
msg = 'xpos and ypos must have the same length'
raise ValueError(msg)
if not (np.all(np.isfinite(xpos)) and np.all(np.isfinite(ypos))):
msg = 'xpos and ypos must contain only finite values'
raise ValueError(msg)
if (xpos.min() < 0 or ypos.min() < 0
or xpos.max() > data.shape[1] - 1
or ypos.max() > data.shape[0] - 1):
msg = 'xpos, ypos values contain points outside the input data'
raise ValueError(msg)
if footprint is None:
if box_size is None:
msg = 'box_size or footprint must be defined'
raise ValueError(msg)
box_size = as_pair('box_size', box_size, lower_bound=(0, 0),
check_odd=True)
footprint = np.ones(box_size, dtype=bool)
else:
footprint = np.asanyarray(footprint, dtype=bool)
if footprint.ndim != 2:
msg = 'footprint must be a 2D array'
raise ValueError(msg)
if not np.any(footprint):
msg = 'footprint must contain at least one True value'
raise ValueError(msg)
if mask is not None and mask.shape != data.shape:
msg = 'mask and data must have the same shape'
raise ValueError(msg)
# Setting error to None is equivalent to no error array, so allow it
# even for centroid functions that do not accept an error keyword
if kwargs.get('error') is None:
kwargs.pop('error', None)
# Allow arbitrary keyword arguments (**kwargs)
spec = inspect.signature(centroid_func)
accepts_var_keyword = any(param.kind == inspect.Parameter.VAR_KEYWORD
for param in spec.parameters.values())
if 'mask' not in spec.parameters and not accepts_var_keyword:
msg = "The input 'centroid_func' must have a 'mask' keyword."
raise ValueError(msg)
if not accepts_var_keyword:
unknown_keys = set(kwargs) - set(spec.parameters)
if unknown_keys:
msg = ('Unrecognized keyword argument(s) for the input '
f"'centroid_func': {sorted(unknown_keys)}")
raise TypeError(msg)
centroid_kwargs = dict(kwargs)
# Save the original error array so that each source independently
# slices the full-image array
error_array = centroid_kwargs.pop('error', None)
if error_array is not None and np.shape(error_array) != data.shape:
msg = 'error and data must have the same shape'
raise ValueError(msg)
# Extract xpeak/ypeak so the original absolute coordinates are
# available for every source. The per-source function below re-adds
# them with the correct cutout offset.
# Remove this block once xpeak and ypeak are fully deprecated.
xpeak_orig = centroid_kwargs.pop('xpeak', None)
ypeak_orig = centroid_kwargs.pop('ypeak', None)
inverted_footprint = np.logical_not(footprint)
def _centroid_source(xypos):
"""
Compute the centroid of the source at the given (x, y)
position.
"""
xp, yp = xypos
slices_large, slices_small = overlap_slices(data.shape,
footprint.shape, (yp, xp))
data_cutout = data[slices_large]
# Trim footprint mask if it has only partial overlap on the data
footprint_mask = inverted_footprint[slices_small]
if mask is not None:
# Combine the input mask cutout and footprint mask
mask_cutout = np.logical_or(mask[slices_large], footprint_mask)
else:
mask_cutout = footprint_mask
if np.all(mask_cutout):
msg = (f'The cutout for the source at ({xp}, {yp}) is completely '
'masked. Please check your input mask and footprint. '
'Also note that footprint must be a small, local '
'footprint.')
raise ValueError(msg)
# Build the per-source keyword arguments from a local copy so
# that no shared state is mutated across sources
src_kwargs = dict(centroid_kwargs)
src_kwargs['mask'] = mask_cutout
if error_array is not None:
src_kwargs['error'] = error_array[slices_large]
# Add xpeak/ypeak with the offset relative to this source's
# cutout.
# Remove this block once xpeak and ypeak are fully deprecated.
if xpeak_orig is not None and ypeak_orig is not None:
src_kwargs['xpeak'] = xpeak_orig - slices_large[1].start
src_kwargs['ypeak'] = ypeak_orig - slices_large[0].start
try:
xcen, ycen = centroid_func(data_cutout, **src_kwargs)
except (ValueError, TypeError) as exc:
msg = f'Centroid failed for source at ({xp}, {yp}): {exc}'
warnings.warn(msg, AstropyUserWarning)
xcen, ycen = np.nan, np.nan
return (xcen + slices_large[1].start,
ycen + slices_large[0].start)
results = [_centroid_source(xypos)
for xypos in zip(xpos, ypos, strict=True)]
results = np.array(results, dtype=float)
return results[:, 0], results[:, 1]
