Merge pull request #2247 from larrybradley/find-peaks-fix

Fix find_peaks edge case with NaN pixels

Author: larrybradley

CHANGES.rst

@@ -198,6 +198,10 @@ Bug Fixes
   - Fixed ``StarFinder`` mutating the input ``kernel`` array during
     normalization. [#2201]
 
+  - Fixed an edge case of ``find_peaks`` not excluding NaN pixel
+    locations from peak detection, which could produce false peaks when
+    the NaN fill value was a local maximum. [#2247]
+
 - ``photutils.isophote``
 
   - ``build_ellipse_model`` now integrates over all angles instead of

photutils/detection/core.py

@@ -49,6 +49,26 @@ class StarFinderBase(metaclass=abc.ABCMeta):
 
     @deprecated_positional_kwargs(since='3.0', until='4.0')
     def __call__(self, data, mask=None):
+        """
+        Find stars in an astronomical image.
+
+        Parameters
+        ----------
+        data : 2D array_like
+            The 2D image array.
+
+        mask : 2D bool array, optional
+            A boolean mask with the same shape as ``data``, where a
+            `True` value indicates the corresponding element of ``data``
+            is masked. Masked pixels are ignored when searching for
+            stars.
+
+        Returns
+        -------
+        table : `~astropy.table.Table` or `None`
+            A table of found stars. If no stars are found then `None` is
+            returned.
+        """
         return self.find_stars(data, mask=mask)
 
     @staticmethod
@@ -60,7 +80,8 @@ def _find_stars(convolved_data, kernel, threshold, *, min_separation=0.0,
         Parameters
         ----------
         convolved_data : 2D array_like
-            The convolved 2D array.
+            The convolved 2D array. Should be NaN-free; any NaN values
+            should be handled before calling this method.
 
         kernel : `_StarFinderKernel` or 2D `~numpy.ndarray`
             The convolution kernel. ``StarFinder`` inputs the kernel
@@ -231,6 +252,7 @@ def __init__(self, data, xypos, kernel, *, n_brightest=None,
         self.xypos = np.atleast_2d(xypos)
         self.n_brightest = n_brightest
         self.peak_max = peak_max
+        self.default_columns = ()
 
         self.id = np.arange(len(self)) + 1
 
@@ -297,7 +319,7 @@ def _get_init_attributes(self):
         This method should be overridden in subclasses.
         """
         return ('data', 'unit', 'kernel', 'n_brightest', 'peak_max',
-                'cutout_shape')
+                'cutout_shape', 'default_columns')
 
     @property
     def _lazyproperties(self):
@@ -583,7 +605,7 @@ def _filter_finite(self, attrs, *, initial_mask=None,
 
         return newcat
 
-    def _filter_bounds(self, bounds, *, peakattr='peak'):
+    def _filter_bounds(self, bounds, *, initial_mask=None, peakattr='peak'):
         """
         Filter the catalog by sharpness, roundness, and peak_max bounds.
 
@@ -595,6 +617,10 @@ def _filter_bounds(self, bounds, *, peakattr='peak'):
             of the form ``(lower_bound, upper_bound)``, or `None` to
             skip filtering for that attribute.
 
+        initial_mask : 1D `~numpy.ndarray` of bool or `None`, optional
+            A pre-existing boolean mask to combine with. If `None`,
+            starts with all `True`.
+
         peakattr : str, optional
             The attribute name for the peak value used for peak_max
             filtering. The default is ``'peak'``.
@@ -604,7 +630,11 @@ def _filter_bounds(self, bounds, *, peakattr='peak'):
         catalog : ``self.__class__`` or `None`
             The filtered catalog, or `None` if no sources remain.
         """
-        mask = np.ones(len(self), dtype=bool)
+        if initial_mask is None:
+            mask = np.ones(len(self), dtype=bool)
+        else:
+            mask = initial_mask.copy()
+
         for attr, range_val in bounds:
             if range_val is None:
                 continue
@@ -670,7 +700,7 @@ def to_table(self, *, columns=None):
 
         table.meta.update(_get_meta())  # keep table.meta type
         if columns is None:
-            if not hasattr(self, 'default_columns'):
+            if not self.default_columns:
                 msg = ('default_columns attribute is not set; either '
                        'pass explicit column names or set '
                        'default_columns in the subclass __init__')

photutils/detection/peakfinder.py

@@ -19,6 +19,71 @@
 __all__ = ['find_peaks']
 
 
+def _verify_ring_candidates(data, peak_mask, needs_verify, footprint_bool,
+                            half, footprint_size):
+    """
+    Verify ring candidates against the exact circular footprint.
+
+    Ring candidates are pixels that are the local maximum within the
+    inscribed box but not in the circumscribed box. These need per-pixel
+    verification against the actual circular footprint.
+
+    Parameters
+    ----------
+    data : 2D `~numpy.ndarray`
+        The 2D image array.
+
+    peak_mask : 2D bool `~numpy.ndarray`
+        Boolean mask to update in place. `True` indicates a confirmed
+        local maximum.
+
+    needs_verify : 2D bool `~numpy.ndarray`
+        Boolean mask of candidate pixels that require verification.
+
+    footprint_bool : 2D bool `~numpy.ndarray`
+        The circular footprint boolean mask.
+
+    half : int
+        Half the footprint size (``footprint_size // 2``), used to
+        center the footprint on each candidate pixel.
+
+    footprint_size : int
+        The size of the circular footprint array along each axis.
+    """
+    y_maybe, x_maybe = needs_verify.nonzero()
+    if len(y_maybe) == 0:
+        return
+
+    ny, nx = data.shape
+    for y, x in zip(y_maybe, x_maybe, strict=True):
+        # Map footprint onto data, clipping to image boundaries
+        y0 = y - half
+        y1 = y0 + footprint_size
+        x0 = x - half
+        x1 = x0 + footprint_size
+
+        dy0, dy1 = max(0, y0), min(ny, y1)
+        dx0, dx1 = max(0, x0), min(nx, x1)
+
+        fy0 = dy0 - y0
+        fy1 = footprint_size - (y1 - dy1)
+        fx0 = dx0 - x0
+        fx1 = footprint_size - (x1 - dx1)
+
+        local = data[dy0:dy1, dx0:dx1]
+        fp_local = footprint_bool[fy0:fy1, fx0:fx1]
+        local_max = local[fp_local].max()
+
+        # Footprint extends beyond image: include cval=0.0
+        if (fy0 > 0 or fy1 < footprint_size or fx0 > 0
+                or fx1 < footprint_size):
+            local_max = max(local_max, 0.0)
+
+        # peak_mask is updated in place
+        if data[y, x] == local_max:
+            peak_mask[y, x] = True
+
+
 def _fast_circular_peaks(data, radius):
     """
     Find pixels that are local maxima within circular regions.
@@ -39,7 +104,9 @@ def _fast_circular_peaks(data, radius):
     Parameters
     ----------
     data : 2D `~numpy.ndarray`
-        The 2D image array (NaN-free).
+        The 2D image array. Must be NaN-free because
+        `~scipy.ndimage.maximum_filter` propagates NaNs, which would
+        corrupt the local-maximum comparisons.
 
     radius : float
         The radius of the circular region in pixels.
@@ -56,7 +123,7 @@ def _fast_circular_peaks(data, radius):
     footprint_size = len(idx)
 
     xx, yy = np.meshgrid(idx, idx)
-    fp_bool = (xx ** 2 + yy ** 2) <= radius_sq
+    footprint_bool = (xx ** 2 + yy ** 2) <= radius_sq
 
     # For even-sized footprints (non-integer radius), scipy's
     # maximum_filter places the center at index ``footprint_size // 2``
@@ -92,36 +159,9 @@ def _fast_circular_peaks(data, radius):
     needs_verify = candidates & ~definite
     peak_mask = definite.copy()
 
-    y_maybe, x_maybe = needs_verify.nonzero()
-    if len(y_maybe) > 0:
-        ny, nx = data.shape
-        for y, x in zip(y_maybe, x_maybe, strict=True):
-
-            # Map footprint onto data, clipping to image boundaries.
-            y0 = y - half
-            y1 = y0 + footprint_size
-            x0 = x - half
-            x1 = x0 + footprint_size
-
-            dy0, dy1 = max(0, y0), min(ny, y1)
-            dx0, dx1 = max(0, x0), min(nx, x1)
-
-            fy0 = dy0 - y0
-            fy1 = footprint_size - (y1 - dy1)
-            fx0 = dx0 - x0
-            fx1 = footprint_size - (x1 - dx1)
-
-            local = data[dy0:dy1, dx0:dx1]
-            fp_local = fp_bool[fy0:fy1, fx0:fx1]
-            local_max = local[fp_local].max()
-
-            # Footprint extends beyond image: include cval=0.0
-            if (fy0 > 0 or fy1 < footprint_size or fx0 > 0
-                    or fx1 < footprint_size):
-                local_max = max(local_max, 0.0)
-
-            if data[y, x] == local_max:
-                peak_mask[y, x] = True
+    # peak_mask is updated in place
+    _verify_ring_candidates(data, peak_mask, needs_verify, footprint_bool,
+                            half, footprint_size)
 
     return peak_mask
 
@@ -262,10 +302,21 @@ def find_peaks(data, threshold, *, box_size=3, footprint=None, mask=None,
 
     A centroiding function can be input via the ``centroid_func``
     keyword to compute centroid coordinates with subpixel precision
-    within the input ``box_size`` or ``footprint``.
+    within the input ``box_size`` or ``footprint``. Note that when
+    ``min_separation`` is used, the centroid region size is determined
+    by ``box_size`` (default 3), not by ``min_separation``.
+
+    The peak detection uses ``mode='constant'`` with ``cval=0.0`` for
+    `~scipy.ndimage.maximum_filter`, which means pixels outside the
+    image boundary are treated as zero. For images with all-negative
+    values, this may suppress legitimate peaks near the borders.
+
+    Any NaN values in the input ``data`` are replaced with the minimum
+    finite value before peak detection, and the corresponding pixels are
+    automatically excluded from the results.
 
-    The output column names (``x_peak``, ``y_peak``,
-    ``peak_value``) differ from the star finder classes (e.g.,
+    The output column names (``x_peak``, ``y_peak``, ``peak_value``)
+    differ from the star finder classes (e.g.,
     `~photutils.detection.DAOStarFinder`), which use ``x_centroid``,
     ``y_centroid``, and ``flux``.
     """
@@ -298,11 +349,14 @@ def find_peaks(data, threshold, *, box_size=3, footprint=None, mask=None,
         border_width = as_pair('border_width', border_width,
                                lower_bound=(0, 1), upper_bound=data.shape)
 
-    # Remove NaN values to avoid runtime warnings
+    # Remove NaN values to avoid runtime warnings and exclude NaN pixels
+    # from peak detection
     nan_mask = np.isnan(data)
     if np.any(nan_mask):
         data = np.copy(data)  # ndarray
         data[nan_mask] = nanmin(data)
+        mask = (nan_mask if mask is None
+                else np.asanyarray(mask) | nan_mask)
 
     # peak_goodmask: good pixels are True
     if min_separation is not None and min_separation > 0:
@@ -318,7 +372,7 @@ def find_peaks(data, threshold, *, box_size=3, footprint=None, mask=None,
 
     # Exclude peaks that are masked
     if mask is not None:
-        mask = np.asanyarray(mask)
+        mask = np.asanyarray(mask, dtype=bool)
         if data.shape != mask.shape:
             msg = 'data and mask must have the same shape'
             raise ValueError(msg)

photutils/detection/tests/test_core.py

@@ -5,6 +5,7 @@
 
 import numpy as np
 import pytest
+from astropy.utils.exceptions import AstropyDeprecationWarning
 
 from photutils.detection import DAOStarFinder
 from photutils.detection.core import (_DEPR_DEFAULT, StarFinderCatalogBase,
@@ -226,7 +227,7 @@ def test_get_init_attributes(self, minimal_catalog_cls):
         xypos = np.array([[5, 5]])
         cat = minimal_catalog_cls(data, xypos, kernel)
         expected = ('data', 'unit', 'kernel', 'n_brightest', 'peak_max',
-                    'cutout_shape')
+                    'cutout_shape', 'default_columns')
         assert cat._get_init_attributes() == expected
 
     def test_lazyproperties_class_cache(self, minimal_catalog_cls):
@@ -520,6 +521,51 @@ def test_getitem_integer_array(self, minimal_catalog_cls):
         assert sub.xypos[0, 0] == 15
         assert sub.xypos[1, 0] == 5
 
+    def test_filter_bounds_none_range(self, minimal_catalog_cls):
+        """
+        Test that _filter_bounds skips filtering when a range is None.
+        """
+        data = np.zeros((21, 21))
+        data[5, 5] = 10.0
+        data[10, 10] = 50.0
+        data[15, 15] = 30.0
+        kernel = np.ones((3, 3))
+        xypos = np.array([[5, 5], [10, 10], [15, 15]])
+        cat = minimal_catalog_cls(data, xypos, kernel)
+        bounds = [('flux', None)]
+        result = cat._filter_bounds(bounds)
+        assert len(result) == 3
+
+    def test_filter_bounds_initial_mask(self, minimal_catalog_cls):
+        """
+        Test _filter_bounds with an initial_mask.
+        """
+        data = np.zeros((21, 21))
+        data[5, 5] = 10.0
+        data[10, 10] = 50.0
+        data[15, 15] = 30.0
+        kernel = np.ones((3, 3))
+        xypos = np.array([[5, 5], [10, 10], [15, 15]])
+        cat = minimal_catalog_cls(data, xypos, kernel)
+        # Pre-exclude the first source
+        initial_mask = np.array([False, True, True])
+        result = cat._filter_bounds([], initial_mask=initial_mask)
+        assert len(result) == 2
+
+    def test_default_columns_preserved_on_slice(self, minimal_catalog_cls):
+        """
+        Test that default_columns is preserved when slicing.
+        """
+        data = np.zeros((21, 21))
+        data[5, 5] = 10.0
+        data[10, 10] = 50.0
+        kernel = np.ones((3, 3))
+        xypos = np.array([[5, 5], [10, 10]])
+        cat = minimal_catalog_cls(data, xypos, kernel)
+        cat.default_columns = ('id', 'x_centroid', 'y_centroid')
+        sub = cat[0]
+        assert sub.default_columns == ('id', 'x_centroid', 'y_centroid')
+
 
 class TestStarFinderBaseCall:
     """
@@ -538,6 +584,18 @@ def test_call_delegates_to_find_stars(self, data):
             np.testing.assert_array_equal(tbl_call[col], tbl_find[col])
 
 
+def test_deprecated_attr(data):
+    """
+    Test that accessing the deprecated attribute on the
+    StarFinderCatalogBase raises an warning.
+    """
+    finder = DAOStarFinder(threshold=5.0, fwhm=2.0)
+    cat = finder._get_raw_catalog(data)
+    match = 'attribute was deprecated'
+    with pytest.warns(AstropyDeprecationWarning, match=match):
+        _ = cat.xcentroid
+
+
 def test_deprecated_default():
     """
     Test repr for _DeprecatedDefault.

photutils/detection/tests/test_daofinder.py

@@ -128,6 +128,20 @@ def test_mask(self, data):
         tbl1 = finder(data, mask=mask)
         assert len(tbl0) > len(tbl1)
 
+    def test_mask_int(self, data):
+        """
+        Test that an integer mask gives the same result as a boolean
+        mask.
+        """
+        finder = DAOStarFinder(threshold=1.0, fwhm=1.5)
+        bool_mask = np.zeros(data.shape, dtype=bool)
+        bool_mask[0:50, :] = True
+        int_mask = bool_mask.astype(int)
+
+        tbl_bool = finder(data, mask=bool_mask)
+        tbl_int = finder(data, mask=int_mask)
+        assert_array_equal(tbl_bool, tbl_int)
+
     def test_xycoords(self, data):
         """
         Test source detection at specified coordinates.