Repository: astropy/photutils

Issue:
New algorithm for robust centroiding
The following is a paper discusses the robustness of algorithms for centroiding stars. Not my area of expertise, but it might be useful to consider if the centroiding algorithm is relevant and if it should be implemented (if it is not already).

[Do fast stellar centroiding methods saturate the Cramér-Rao lower bound?](https://arxiv.org/abs/1610.05873 )
Mohammadjavad Vakili, David W. Hogg
(Submitted on 19 Oct 2016)

    One of the most demanding tasks in astronomical image processing---in terms of precision---is the centroiding of stars. Upcoming large surveys are going to take images of billions of point sources, including many faint stars, with short exposure times. Real-time estimation of the centroids of stars is crucial for real-time PSF estimation, and maximal precision is required for measurements of proper motion.
    The fundamental Cramér-Rao lower bound sets a limit on the root-mean-squared-error achievable by optimal estimators. In this work, we aim to compare the performance of various centroiding methods, in terms of saturating the bound, when they are applied to relatively low signal-to-noise ratio unsaturated stars assuming zero-mean constant Gaussian noise. In order to make this comparison, we present the ratio of the root-mean-squared-errors of these estimators to their corresponding Cramér-Rao bound as a function of the signal-to-noise ratio and the full-width at half-maximum of faint stars.
    We discuss two general circumstances in centroiding of faint stars: (i) when we have a good estimate of the PSF, (ii) when we do not know the PSF. In the case that we know the PSF, we show that a fast polynomial centroiding after smoothing the image by the PSF can be as efficient as the maximum-likelihood estimator at saturating the bound. In the case that we do not know the PSF, we demonstrate that although polynomial centroiding is not as optimal as PSF profile fitting, it comes very close to saturating the Cramér-Rao lower bound in a wide range of conditions. We also show that the moment-based method of center-of-light never comes close to saturating the bound, and thus it does not deliver reliable estimates of centroids. 

Discussion:

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index dbca1405..9b69b3fa 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -19,6 +19,11 @@ New Features
   - Added ``coverage_mask`` and ``fill_value`` keyword options to
     ``Background2D``. [#1061]
 
+- ``photutils.centroids``
+
+  - Added quadratic centroid estimator function
+    (``centroid_quadratic``). [#1067]
+
 Bug Fixes
 ^^^^^^^^^
 
diff --git a/docs/centroids.rst b/docs/centroids.rst
index 600692cd..e6100a97 100644
--- a/docs/centroids.rst
+++ b/docs/centroids.rst
@@ -10,6 +10,9 @@ centroid of a single source:
 * :func:`~photutils.centroids.centroid_com`: Calculates the object
   "center of mass" from 2D image moments.
 
+* :func:`~photutils.centroids.centroid_quadratic`: Calculates the
+  centroid by fitting a 2D quadratic polynomial to the data.
+
 * :func:`~photutils.centroids.centroid_1dg`: Calculates the centroid
   by fitting 1D Gaussians to the marginal ``x`` and ``y``
   distributions of the data.
@@ -18,13 +21,13 @@ centroid of a single source:
   by fitting a 2D Gaussian to the 2D distribution of the data.
 
 Masks can be input into each of these functions to mask bad pixels.
-Error arrays can be input into the two fitting methods to weight the
-fits.
+Error arrays can be input into the two Gaussian fitting methods to
+weight the fits.
 
 To calculate the centroids of many sources in an image, use the
 :func:`~photutils.centroids.centroid_sources` function. This function
 can be used with any of the above centroiding functions or a custom
-user-defined centroiding function.  See zzzzz
+user-defined centroiding function.
 
 
 Getting Started
diff --git a/photutils/centroids/core.py b/photutils/centroids/core.py
index 5c069d79..8e4afd79 100644
--- a/photutils/centroids/core.py
+++ b/photutils/centroids/core.py
@@ -10,7 +10,10 @@ from astropy.nddata.utils import overlap_slices
 from astropy.utils.exceptions import AstropyUserWarning
 import numpy as np
 
-__all__ = ['centroid_com', 'centroid_sources', 'centroid_epsf']
+from ..utils._round import _py2intround
+
+__all__ = ['centroid_com', 'centroid_quadratic', 'centroid_sources',
+           'centroid_epsf']
 
 
 def centroid_com(data, mask=None, oversampling=1):
@@ -72,6 +75,181 @@ def centroid_com(data, mask=None, oversampling=1):
                      for axis in range(data.ndim)])[::-1]
 
 
+def centroid_quadratic(data, xpeak=None, ypeak=None, fit_boxsize=5,
+                       search_boxsize=None, mask=None):
+    """
+    Calculate the centroid of an n-dimensional array by fitting a 2D
+    quadratic polynomial.
+
+    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 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 ``xmax`` and ``ymax`` 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``.
+
+    Parameters
+    ----------
+    image_data : numpy.ndarray
+        Image data.
+
+    xpeak, ypeak : float or `None`, optional
+        The initial guess of the position of the centroid. When both
+        ``xpeak`` and ``ypeak`` are `None` then the position of the
+        maximum value in the input ``data`` will be used as the initial
+        guess.
+
+    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 should be in
+        ``(ny, nx)`` order. If ``fit_boxsize`` is a scalar then a square
+        box of size ``fit_boxsize`` will be used.
+
+    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 `None`. If
+        ``fit_boxsize`` has two elements, they should be in ``(ny, nx)``
+        order. If ``fit_boxsize`` is a scalar then a square box of size
+        ``fit_boxsize`` will be used.  This parameter is ignored when
+        ``xmax`` and ``ymax`` are both `None`.  In that case, the entire
+        array is search for the maximum value.
+
+    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.
+        Masked data are excluded from calculations.
+
+    Returns
+    -------
+    centroid : `~numpy.ndarray`
+        The ``x, y`` coordinates of the centroid.
+    """
+    if ((xpeak is None and ypeak is not None)
+            or (xpeak is not None and ypeak is None)):
+        raise ValueError('xpeak and ypeak must both be input or "None"')
+
+    data = np.asanyarray(data, dtype=float)
+    ny, nx = data.shape
+
+    badmask = ~np.isfinite(data)
+    if np.any(badmask):
+        warnings.warn('Input data contains non-finite values (e.g., NaN or '
+                      'inf) that were automatically masked.',
+                      AstropyUserWarning)
+        data[badmask] = np.nan
+
+    if mask is not None:
+        if data.shape != mask.shape:
+            raise ValueError('data and mask must have the same shape.')
+        data[mask] = np.nan
+
+    fit_boxsize = _process_boxsize(fit_boxsize, data.shape)
+
+    if np.product(fit_boxsize) < 6:
+        raise ValueError('fit_boxsize is too small.  6 values are required '
+                         'to fit a 2D quadratic polynomial.')
+
+    if xpeak is None:  # and ypeak too
+        yidx, xidx = np.unravel_index(np.nanargmax(data), data.shape)
+    else:
+        xidx = _py2intround(xpeak)
+        yidx = _py2intround(ypeak)
+
+        if search_boxsize is not None:
+            search_boxsize = _process_boxsize(search_boxsize, data.shape)
+
+            slc_data, slc_cutout = 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
+
+    # if peak is at the edge of the data, return the position of the maximum
+    if xidx == 0 or xidx == nx - 1 or yidx == 0 or yidx == ny - 1:
+        warnings.warn('maximum value is at the edge of the data and its '
+                      'position was returned; no quadratic fit was '
+                      'performed', AstropyUserWarning)
+        return np.array((xidx, yidx), dtype=float)
+
+    # extract the fitting region
+    slc_data, slc_cutout = 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:
+        warnings.warn('at least 6 unmasked data points are required to '
+                      'perform a 2D quadratic fit',
+                      AstropyUserWarning)
+        return np.array((np.nan, np.nan))
+
+    # fit a 2D quadratic polynomial to the fitting region
+    xi = np.arange(xidx0, xidx1)
+    yi = np.arange(yidx0, yidx1)
+    x, y = np.meshgrid(xi, yi)
+    x = x.ravel()
+    y = y.ravel()
+    coeff_matrix = np.vstack((np.ones_like(x), x, y, x * y, x * x, y * y)).T
+
+    try:
+        c = np.linalg.lstsq(coeff_matrix, cutout, rcond=None)[0]
+    except np.linalg.LinAlgError:
+        warnings.warn('quadratic fit failed', 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 det <= 0 or ((c20 > 0.0 and c02 >= 0.0) or (c20 >= 0.0 and c02 > 0.0)):
+        warnings.warn('quadratic fit does not have a maximum',
+                      AstropyUserWarning)
+        return np.array((np.nan, np.nan))
+
+    xm = (c01 * c11 - 2.0 * c02 * c10) / det
+    ym = (c10 * c11 - 2.0 * c20 * c01) / det
+    if xm > 0.0 and xm < (nx - 1.0) and ym > 0.0 and ym < (ny - 1.0):
+        xycen = np.array((xm, ym), dtype=float)
+    else:
+        warnings.warn('quadratic polynomial maximum value falls outside '
+                      'of the image', AstropyUserWarning)
+        return np.array((np.nan, np.nan))
+
+    return xycen
+
+
+def _process_boxsize(box_size, data_shape):
+    box_size = np.round(np.atleast_1d(box_size)).astype(int)
+    if len(box_size) == 1:
+        box_size = np.repeat(box_size, 2)
+    if len(box_size) > 2:
+        raise ValueError('box size must contain only 1 or 2 values')
+    if np.any(box_size < 0):
+        raise ValueError('box size must be >= 0')
+    # box_size cannot be larger than the data shape
+    box_size = (min(box_size[0], data_shape[0]),
+                min(box_size[1], data_shape[1]))
+    return box_size
+
+
 def centroid_sources(data, xpos, ypos, box_size=11, footprint=None,
                      error=None, mask=None, centroid_func=centroid_com):
     """
diff --git a/photutils/psf/epsf.py b/photutils/psf/epsf.py
index d9c9624e..8c1491eb 100644
--- a/photutils/psf/epsf.py
+++ b/photutils/psf/epsf.py
@@ -19,6 +19,7 @@ import numpy as np
 from .epsf_stars import EPSFStar, EPSFStars, LinkedEPSFStar
 from .models import EPSFModel
 from ..centroids import centroid_epsf
+from ..utils._round import _py2intround
 
 try:
     import bottleneck  # pylint: disable=W0611
@@ -874,23 +875,6 @@ class EPSFBuilder:
         return epsf, stars
 
 
-def _py2intround(a):
-    """
-    Round the input to the nearest integer.
-
-    If two integers are equally close, rounding is done away from 0.
-    """
-
-    data = np.asanyarray(a)
-    value = np.where(data >= 0, np.floor(data + 0.5),
-                     np.ceil(data - 0.5)).astype(int)
-
-    if not hasattr(a, '__iter__'):
-        value = np.asscalar(value)
-
-    return value
-
-
 def _interpolate_missing_data(data, mask, method='cubic'):
     """
     Interpolate missing data as identified by the ``mask`` keyword.
diff --git a/photutils/utils/_round.py b/photutils/utils/_round.py
new file mode 100644
index 00000000..31daea71
--- /dev/null
+++ b/photutils/utils/_round.py
@@ -0,0 +1,23 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""
+This module provides tools to round numpy arrays.
+"""
+
+import numpy as np
+
+
+def _py2intround(a):
+    """
+    Round the input to the nearest integer.
+
+    If two integers are equally close, rounding is done away from 0.
+    """
+
+    data = np.asanyarray(a)
+    value = np.where(data >= 0, np.floor(data + 0.5),
+                     np.ceil(data - 0.5)).astype(int)
+
+    if not hasattr(a, '__iter__'):
+        value = np.asscalar(value)
+
+    return value


Test patch:
diff --git a/photutils/centroids/tests/test_core.py b/photutils/centroids/tests/test_core.py
index 0761b6c5..9dd01e83 100644
--- a/photutils/centroids/tests/test_core.py
+++ b/photutils/centroids/tests/test_core.py
@@ -12,7 +12,8 @@ import numpy as np
 from numpy.testing import assert_allclose
 import pytest
 
-from ..core import centroid_com, centroid_sources, centroid_epsf
+from ..core import (centroid_com, centroid_quadratic, centroid_sources,
+                    centroid_epsf)
 from ..gaussian import centroid_1dg
 from ...psf import IntegratedGaussianPRF
 
@@ -47,6 +48,9 @@ def test_centroid_com(x_std, y_std, theta):
     xc, yc = centroid_com(data)
     assert_allclose((xc, yc), (XCEN, YCEN), rtol=0, atol=1.e-3)
 
+    xc, yc = centroid_quadratic(data)
+    assert_allclose((xc, yc), (XCEN, YCEN), rtol=0, atol=0.015)
+
     # test with mask
     mask = np.zeros(data.shape, dtype=bool)
     data[10, 10] = 1.e5
@@ -54,6 +58,9 @@ def test_centroid_com(x_std, y_std, theta):
     xc, yc = centroid_com(data, mask=mask)
     assert_allclose((xc, yc), (XCEN, YCEN), rtol=0, atol=1.e-3)
 
+    xc, yc = centroid_quadratic(data, mask=mask)
+    assert_allclose((xc, yc), (XCEN, YCEN), rtol=0, atol=0.015)
+
     # test with oversampling
     for oversampling in [4, (4, 6)]:
         if not hasattr(oversampling, '__len__'):
@@ -91,6 +98,13 @@ def test_centroid_com_nan_withmask(use_mask):
     if nwarn == 1:
         assert issubclass(warnlist[0].category, AstropyUserWarning)
 
+    with warnings.catch_warnings(record=True) as warnlist:
+        xc, yc = centroid_quadratic(data, mask=mask)
+        assert_allclose(xc, xc_ref, rtol=0, atol=0.15)
+    assert len(warnlist) == 1  # always warns because of NaN
+    if nwarn == 1:
+        assert issubclass(warnlist[0].category, AstropyUserWarning)
+
 
 @pytest.mark.skipif('not HAS_SCIPY')
 def test_centroid_com_invalid_inputs():
@@ -102,10 +116,73 @@ def test_centroid_com_invalid_inputs():
         centroid_com(data, oversampling=-1)
 
 
+@pytest.mark.skipif('not HAS_SCIPY')
+def test_centroid_quadratic_xypeak():
+    data = np.zeros((11, 11))
+    data[4:7, 4:7] = 10
+    data[5, 5] = 100
+    xycen1 = centroid_quadratic(data)
+    xycen2 = centroid_quadratic(data, xpeak=5, ypeak=5)
+    xycen3 = centroid_quadratic(data, xpeak=5, ypeak=5, search_boxsize=3)
+    assert_allclose(xycen1, xycen2)
+    assert_allclose(xycen1, xycen3)
+
+
+@pytest.mark.skipif('not HAS_SCIPY')
+def test_centroid_quadratic_npts():
+    data = np.zeros((3, 3))
+    data[1, 1] = 1
+    mask = np.zeros(data.shape, dtype=bool)
+    mask[0, :] = True
+    mask[2, :] = True
+    with warnings.catch_warnings(record=True) as warnlist:
+        centroid_quadratic(data, mask=mask)
+    assert issubclass(warnlist[0].category, AstropyUserWarning)
+
+
+@pytest.mark.skipif('not HAS_SCIPY')
+def test_centroid_quadratic_invalid_inputs():
+    data = np.zeros((4, 4))
+    mask = np.zeros((2, 2), dtype=bool)
+    with pytest.raises(ValueError):
+        centroid_quadratic(data, xpeak=3, ypeak=None)
+    with pytest.raises(ValueError):
+        centroid_quadratic(data, xpeak=None, ypeak=3)
+    with pytest.raises(ValueError):
+        centroid_quadratic(data, fit_boxsize=(2, 2, 2))
+    with pytest.raises(ValueError):
+        centroid_quadratic(data, fit_boxsize=(-2, 2))
+    with pytest.raises(ValueError):
+        centroid_quadratic(data, fit_boxsize=(2, 2))
+    with pytest.raises(ValueError):
+        centroid_quadratic(data, mask=mask)
+
+
+@pytest.mark.skipif('not HAS_SCIPY')
+def test_centroid_quadratic_edge():
+    data = np.zeros((11, 11))
+    data[1, 1] = 100
+    data[9, 9] = 100
+
+    xycen = centroid_quadratic(data, xpeak=1, ypeak=1, fit_boxsize=5)
+    assert_allclose(xycen, (0.923077, 0.923077))
+
+    xycen = centroid_quadratic(data, xpeak=9, ypeak=9, fit_boxsize=5)
+    assert_allclose(xycen, (9.076923, 9.076923))
+
+    data = np.zeros((5, 5))
+    data[0, 0] = 100
+    with warnings.catch_warnings(record=True) as warnlist:
+        xycen = centroid_quadratic(data)
+    assert_allclose(xycen, (0, 0))
+    assert issubclass(warnlist[0].category, AstropyUserWarning)
+
+
 @pytest.mark.skipif('not HAS_SCIPY')
 def test_centroid_sources():
+    theta = np.pi / 6.
     model = Gaussian2D(2.4, XCEN, YCEN, x_stddev=3.2, y_stddev=5.7,
-                       theta=np.pi/6.)
+                       theta=theta)
     y, x = np.mgrid[0:50, 0:47]
     data = model(x, y)
     error = np.ones(data.shape, dtype=float)