Merge pull request #1943 from larrybradley/fix-finder-border

Fix issues with exclude_border option in star finders

CHANGES.rst

@@ -18,6 +18,14 @@ Bug Fixes
     using the ``BkgIDWInterpolator`` and when any mesh was excluded,
     e.g., due to an input mask. [#1940]
 
+- ``photutils.detection``
+
+  - Fixed a bug in the star finders (``DAOStarFinder``,
+    ``IRAFStarFinder``, and ``StarFinder``) when
+    ``exclude_border=True``. Also, fixed an issue with
+    ``exclude_border=True`` where if all sources were in the border
+    region then an error would be raised. [#1943]
+
 
 2.0.1 (2024-10-16)
 ------------------

photutils/detection/core.py

@@ -39,8 +39,9 @@ def _find_stars(convolved_data, kernel, threshold, *, min_separation=0.0,
         convolved_data : 2D array_like
             The convolved 2D array.
 
-        kernel : `_StarFinderKernel`
-            The convolution kernel.
+        kernel : `_StarFinderKernel` or 2D `~numpy.ndarray`
+            The convolution kernel. ``StarFinder`` inputs the kernel
+            as a 2D array.
 
         threshold : float
             The absolute image value above which to select sources. This
@@ -83,46 +84,29 @@ class multiplied by the kernel relerr. If ``convolved_data``
             footprint = np.array((xx**2 + yy**2) <= min_separation**2,
                                  dtype=int)
 
-        # pad the convolved data and mask by half the kernel size (or
-        # x/y radius) to allow for detections near the edges
-        if isinstance(kernel, np.ndarray):
-            ypad = (kernel.shape[0] - 1) // 2
-            xpad = (kernel.shape[1] - 1) // 2
+        # define the border exclusion region
+        if exclude_border:
+            if isinstance(kernel, np.ndarray):
+                yborder = (kernel.shape[0] - 1) // 2
+                xborder = (kernel.shape[1] - 1) // 2
+            else:
+                yborder = kernel.yradius
+                xborder = kernel.xradius
+            border_width = max(xborder, yborder)
         else:
-            ypad = kernel.yradius
-            xpad = kernel.xradius
-
-        if not exclude_border:
-            pad = ((ypad, ypad), (xpad, xpad))
-            pad_mode = 'constant'
-            convolved_data = np.pad(convolved_data, pad, mode=pad_mode,
-                                    constant_values=0.0)
-            if mask is not None:
-                mask = np.pad(mask, pad, mode=pad_mode, constant_values=False)
+            border_width = 0
 
         # find local peaks in the convolved data
         # suppress any NoDetectionsWarning from find_peaks
         with warnings.catch_warnings():
             warnings.filterwarnings('ignore', category=NoDetectionsWarning)
             tbl = find_peaks(convolved_data, threshold, footprint=footprint,
-                             mask=mask)
+                             mask=mask, border_width=border_width)
 
         if tbl is None:
             return None
 
-        if exclude_border:
-            xmax = convolved_data.shape[1] - xpad
-            ymax = convolved_data.shape[0] - ypad
-            mask = ((tbl['x_peak'] > xpad) & (tbl['y_peak'] > ypad)
-                    & (tbl['x_peak'] < xmax) & (tbl['y_peak'] < ymax))
-            tbl = tbl[mask]
-
-        xpos, ypos = tbl['x_peak'], tbl['y_peak']
-        if not exclude_border:
-            xpos -= xpad
-            ypos -= ypad
-
-        return np.transpose((xpos, ypos))
+        return np.transpose((tbl['x_peak'], tbl['y_peak']))
 
     @abc.abstractmethod
     def find_stars(self, data, mask=None):

photutils/detection/irafstarfinder.py

@@ -416,7 +416,7 @@ def sky(self):
         skymask = ~self.kernel.mask.astype(bool)  # 1=sky, 0=obj
         nsky = np.count_nonzero(skymask)
         axis = (1, 2)
-        if nsky == 0.0:
+        if nsky == 0.0:  # pragma: no cover
             sky = (np.max(self.cutout_data_nosub, axis=axis)
                    - np.max(self.cutout_convdata, axis=axis))
         else:
@@ -452,7 +452,7 @@ def cutout_data(self):
         return data
 
     @lazyproperty
-    def cutout_convdata(self):
+    def cutout_convdata(self):  # pragma: no cover
         return self.make_cutouts(self.convolved_data)
 
     @lazyproperty

photutils/detection/peakfinder.py

@@ -19,7 +19,7 @@
 
 
 def find_peaks(data, threshold, *, box_size=3, footprint=None, mask=None,
-               border_width=None, npeaks=np.inf, centroid_func=None,
+               border_width=0, npeaks=np.inf, centroid_func=None,
                error=None, wcs=None):
     """
     Find local peaks in an image that are above a specified threshold
@@ -74,9 +74,9 @@ def find_peaks(data, threshold, *, box_size=3, footprint=None, mask=None,
         A boolean mask with the same shape as ``data``, where a `True`
         value indicates the corresponding element of ``data`` is masked.
 
-    border_width : bool, optional
+    border_width : int, optional
         The width in pixels to exclude around the border of the
-        ``data``.
+        ``data``. Must be an non-negative integer.
 
     npeaks : int, optional
         The maximum number of peaks to return. When the number of
@@ -126,7 +126,6 @@ def find_peaks(data, threshold, *, box_size=3, footprint=None, mask=None,
     arrays, unit = process_quantities((data, threshold, error),
                                       ('data', 'threshold', 'error'))
     data, threshold, error = arrays
-
     data = np.asanyarray(data)
 
     if np.all(data == data.flat[0]):
@@ -140,6 +139,11 @@ def find_peaks(data, threshold, *, box_size=3, footprint=None, mask=None,
             raise ValueError('A threshold array must have the same shape as '
                              'the input data.')
 
+    if border_width < 0:
+        raise ValueError('border_width must be a non-negative integer.')
+    if int(border_width) != border_width:
+        raise ValueError('border_width must be an integer.')
+
     # remove NaN values to avoid runtime warnings
     nan_mask = np.isnan(data)
     if np.any(nan_mask):
@@ -161,7 +165,7 @@ def find_peaks(data, threshold, *, box_size=3, footprint=None, mask=None,
             raise ValueError('data and mask must have the same shape')
         peak_goodmask = np.logical_and(peak_goodmask, ~mask)
 
-    if border_width is not None:
+    if border_width > 0:
         for i in range(peak_goodmask.ndim):
             peak_goodmask = peak_goodmask.swapaxes(0, i)
             peak_goodmask[:border_width] = False

photutils/detection/starfinder.py

@@ -96,7 +96,8 @@ def _get_raw_catalog(self, data, *, mask=None):
         kernel = self.kernel
         kernel /= np.max(kernel)  # normalize max value to 1.0
         denom = np.sum(kernel**2) - (np.sum(kernel)**2 / kernel.size)
-        kernel = (kernel - np.sum(kernel) / kernel.size) / denom
+        if denom > 0:
+            kernel = (kernel - np.sum(kernel) / kernel.size) / denom
 
         convolved_data = _filter_data(data, kernel, mode='constant',
                                       fill_value=0.0,

photutils/detection/tests/test_irafstarfinder.py

@@ -9,6 +9,7 @@
 from numpy.testing import assert_array_equal
 
 from photutils.detection import IRAFStarFinder
+from photutils.psf import CircularGaussianPRF
 from photutils.utils.exceptions import NoDetectionsWarning
 
 
@@ -171,3 +172,30 @@ def test_single_detected_source(self, data):
         assert cat.isscalar
         flux = cat.flux[0]  # evaluate the flux so it can be sliced
         assert cat[0].flux == flux
+
+    def test_all_border_sources(self):
+        model1 = CircularGaussianPRF(flux=100, x_0=1, y_0=1, fwhm=2)
+        model2 = CircularGaussianPRF(flux=100, x_0=50, y_0=50, fwhm=2)
+        model3 = CircularGaussianPRF(flux=100, x_0=30, y_0=30, fwhm=2)
+
+        threshold = 1
+        yy, xx = np.mgrid[:51, :51]
+        data = model1(xx, yy)
+
+        # test single source within the border region
+        finder = IRAFStarFinder(threshold=threshold, fwhm=2.0, roundlo=-0.1,
+                                exclude_border=True)
+        with pytest.warns(NoDetectionsWarning):
+            tbl = finder(data)
+        assert tbl is None
+
+        # test multiple sources all within the border region
+        data += model2(xx, yy)
+        with pytest.warns(NoDetectionsWarning):
+            tbl = finder(data)
+        assert tbl is None
+
+        # test multiple sources with some within the border region
+        data += model3(xx, yy)
+        tbl = finder(data)
+        assert len(tbl) == 1

photutils/detection/tests/test_peakfinder.py

@@ -91,6 +91,13 @@ def test_border_width(self, data):
         tbl1 = find_peaks(data, 0.1, box_size=3, border_width=25)
         assert len(tbl1) < len(tbl0)
 
+        match = 'border_width must be a non-negative integer'
+        with pytest.raises(ValueError, match=match):
+            find_peaks(data, 0.1, box_size=3, border_width=-1)
+        match = 'border_width must be an integer'
+        with pytest.raises(ValueError, match=match):
+            find_peaks(data, 0.1, box_size=3, border_width=3.1)
+
     def test_box_size_int(self, data):
         """
         Test non-integer box_size.

photutils/detection/tests/test_starfinder.py

@@ -46,6 +46,18 @@ def test_inputs(self, kernel):
         with pytest.raises(ValueError, match=match):
             StarFinder(1, kernel, brightest=3.1)
 
+    def test_exclude_border(self, data, kernel):
+        data = np.zeros((12, 12))
+        data[0:2, 0:2] = 1
+        data[9:12, 9:12] = 1
+        kernel = np.ones((3, 3))
+
+        finder0 = StarFinder(1, kernel, exclude_border=False)
+        finder1 = StarFinder(1, kernel, exclude_border=True)
+        tbl0 = finder0(data)
+        tbl1 = finder1(data)
+        assert len(tbl0) > len(tbl1)
+
     def test_nosources(self, data, kernel):
         match = 'No sources were found'
         with pytest.warns(NoDetectionsWarning, match=match):