# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Tools for centroiding sources using Gaussians.
"""
import warnings
import numpy as np
from astropy.modeling.fitting import TRFLSQFitter
from astropy.modeling.models import Gaussian1D, Gaussian2D
from astropy.utils.exceptions import AstropyUserWarning
from photutils.centroids._utils import (_gaussian1d_moments,
_gaussian2d_moments,
_validate_gaussian_inputs)
from photutils.utils._deprecation import deprecated_positional_kwargs
from photutils.utils._quantity_helpers import process_quantities
__all__ = ['centroid_1dg', 'centroid_2dg']
[docs]
@deprecated_positional_kwargs(since='3.0', until='4.0')
def centroid_1dg(data, error=None, mask=None):
"""
Calculate the centroid of a 2D array by fitting 1D Gaussians to the
marginal ``x`` and ``y`` distributions of the array.
Non-finite values (e.g., NaN or inf) in the ``data`` or ``error``
arrays 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`.
Masked pixels are excluded by substituting zero into the
marginal sums, and the fit weights are zeroed only for
fully-masked rows or columns. A partially-masked row or column
near the source peak therefore distorts the corresponding
marginal distribution and can bias the centroid. Consider using
`~photutils.centroids.centroid_2dg`, which excludes individual
masked pixels from the fit, when isolated masked or non-finite
pixels fall near the source peak.
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.
error : 2D `~numpy.ndarray`, optional
The 2D array of the 1-sigma errors of the input ``data``.
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``.
Returns
-------
centroid : `~numpy.ndarray`
The ``x, y`` coordinates of the centroid.
Examples
--------
>>> import numpy as np
>>> from photutils.datasets import make_4gaussians_image
>>> from photutils.centroids import centroid_1dg
>>> data = make_4gaussians_image()
>>> data -= np.median(data[0:30, 0:125])
>>> data = data[40:80, 70:110]
>>> x1, y1 = centroid_1dg(data)
>>> print(np.array((x1, y1)))
[19.96553246 20.04952841]
.. plot::
import matplotlib.pyplot as plt
import numpy as np
from photutils.centroids import centroid_1dg
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_1dg(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, error), _ = process_quantities((data, error), ('data', 'error'))
data, mask, error = _validate_gaussian_inputs(data, mask, error)
if error is not None:
error_squared = error**2
xy_error = [np.sqrt(np.sum(error_squared, axis=i)) for i in (0, 1)]
xy_weights = [1.0 / xy_err.clip(min=1.0e-30) for xy_err in xy_error]
else:
xy_weights = [np.ones(data.shape[i]) for i in (1, 0)]
# Assign zero weight where an entire row or column is masked
if np.any(mask):
bad_idx = [np.all(mask, axis=i) for i in (0, 1)]
for i in (0, 1):
xy_weights[i][bad_idx[i]] = 0.0
xy_data = [np.sum(data, axis=i) for i in (0, 1)]
# A zero-sum or constant marginal distribution has undefined
# moment estimates (0/0) and cannot be fit with a 1D Gaussian.
for data_i in xy_data:
if np.sum(data_i) == 0 or np.ptp(data_i) == 0:
msg = ('Input data must have a nonzero sum and non-constant '
'values to fit 1D Gaussians.')
raise ValueError(msg)
# Gaussian1D stddev is bounded to be strictly positive
fitter = TRFLSQFitter()
centroid = []
for (data_i, weights_i) in zip(xy_data, xy_weights, strict=True):
params_init = _gaussian1d_moments(data_i)
g_init = Gaussian1D(*params_init)
x = np.arange(data_i.size)
g_fit = fitter(g_init, x, data_i, weights=weights_i)
centroid.append(g_fit.mean.value)
return np.array(centroid)
[docs]
@deprecated_positional_kwargs(since='3.0', until='4.0')
def centroid_2dg(data, error=None, mask=None):
"""
Calculate the centroid of a 2D array by fitting a 2D Gaussian to the
array.
Non-finite values (e.g., NaN or inf) in the ``data`` or ``error``
arrays 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 : 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.
error : 2D `~numpy.ndarray`, optional
The 2D array of the 1-sigma errors of the input ``data``.
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``.
Returns
-------
centroid : `~numpy.ndarray`
The ``x, y`` coordinates of the centroid.
Examples
--------
>>> import numpy as np
>>> from photutils.datasets import make_4gaussians_image
>>> from photutils.centroids import centroid_2dg
>>> data = make_4gaussians_image()
>>> data -= np.median(data[0:30, 0:125])
>>> data = data[40:80, 70:110]
>>> x1, y1 = centroid_2dg(data)
>>> print(np.array((x1, y1)))
[19.9851944 20.01490157]
.. plot::
import matplotlib.pyplot as plt
import numpy as np
from photutils.centroids import centroid_2dg
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_2dg(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, error), _ = process_quantities((data, error), ('data', 'error'))
data, mask, error = _validate_gaussian_inputs(data, mask, error)
if np.count_nonzero(~mask) < 6:
msg = ('Input data must have at least 6 unmasked values to fit a '
'2D Gaussian.')
raise ValueError(msg)
# Subtract the minimum of the data to make the data values positive.
# Moments from negative data values can yield undefined Gaussian
# parameters, e.g., x_stddev and y_stddev.
shifted = data - np.min(data)
if np.sum(shifted) == 0:
msg = ('Input data must have non-constant values to fit a '
'2D Gaussian.')
raise ValueError(msg)
if error is not None:
weights = 1.0 / error.clip(min=1.0e-30)
else:
weights = np.ones(data.shape)
# Assign zero weight to masked pixels
if np.any(mask):
weights[mask] = 0.0
amplitude, x_mean, y_mean, x_stddev, y_stddev, theta = _gaussian2d_moments(
shifted)
g_init = Gaussian2D(amplitude=amplitude,
x_mean=x_mean,
y_mean=y_mean,
x_stddev=x_stddev,
y_stddev=y_stddev,
theta=theta)
fitter = TRFLSQFitter()
y, x = np.indices(data.shape)
gfit = fitter(g_init, x, y, data, weights=weights)
# TRFLSQFitter stores the scipy least_squares result object in
# fit_info. Success is False when the optimizer terminated without
# satisfying a convergence criterion (e.g., the maximum number of
# function evaluations was exceeded). Inspecting fit_info instead
# of capturing warnings avoids mutating the process-global warnings
# state, which is not thread-safe.
if not fitter.fit_info.success:
msg = 'The fit may not have converged. Please check your results.'
warnings.warn(msg, AstropyUserWarning)
return np.array([gfit.x_mean.value, gfit.y_mean.value])
