Merge pull request #1663 from larrybradley/daofind-minsep

Add min_separation keyword to DAOStarFinder and IRAFStarFinder

CHANGES.rst

@@ -9,6 +9,11 @@ General
 New Features
 ^^^^^^^^^^^^
 
+- ``photutils.detection``
+
+  - Added a ``min_separation`` keyword to ``DAOStarFinder`` and
+    ``IRAFStarFinder``. [#1663]
+
 - ``photutils.morphology``
 
   - Added a ``wcs`` keyword to ``data_properties``. [#1648]

photutils/detection/core.py

@@ -69,7 +69,7 @@ def _find_stars(convolved_data, kernel, threshold, *, min_separation=0.0,
         .. _starfind: https://iraf.net/irafhelp.php?val=starfind
         """
         # define a local footprint for the peak finder
-        if min_separation == 0:  # daofind
+        if min_separation == 0.0:  # DAOStarFinder
             if isinstance(kernel, np.ndarray):
                 footprint = np.ones(kernel.shape)
             else:

photutils/detection/daofinder.py

@@ -119,11 +119,15 @@ class DAOStarFinder(StarFinderBase):
             pixel values are negative. Therefore, setting ``peakmax`` to a
             non-positive value would result in exclusion of all objects.
 
-    xycoords : `None` or Nx2 `~numpy.ndarray`
+    xycoords : `None` or Nx2 `~numpy.ndarray`, optional
         The (x, y) pixel coordinates of the approximate centroid
         positions of identified sources. If ``xycoords`` are input, the
         algorithm will skip the source-finding step.
 
+    min_separation : float, optional
+        The minimum separation (in pixels) for detected objects. Note
+        that large values may result in long run times.
+
     See Also
     --------
     IRAFStarFinder
@@ -157,7 +161,8 @@ class DAOStarFinder(StarFinderBase):
     def __init__(self, threshold, fwhm, ratio=1.0, theta=0.0,
                  sigma_radius=1.5, sharplo=0.2, sharphi=1.0, roundlo=-1.0,
                  roundhi=1.0, sky=0.0, exclude_border=False,
-                 brightest=None, peakmax=None, xycoords=None):
+                 brightest=None, peakmax=None, xycoords=None,
+                 min_separation=0.0):
 
         if not np.isscalar(threshold):
             raise TypeError('threshold must be a scalar value.')
@@ -179,6 +184,10 @@ def __init__(self, threshold, fwhm, ratio=1.0, theta=0.0,
         self.brightest = self._validate_brightest(brightest)
         self.peakmax = peakmax
 
+        if min_separation < 0:
+            raise ValueError('min_separation must be >= 0')
+        self.min_separation = min_separation
+
         if xycoords is not None:
             xycoords = np.asarray(xycoords)
             if xycoords.ndim != 2 or xycoords.shape[1] != 2:
@@ -208,6 +217,7 @@ def _get_raw_catalog(self, data, mask=None):
         if self.xycoords is None:
             xypos = self._find_stars(convolved_data, self.kernel,
                                      self.threshold_eff, mask=mask,
+                                     min_separation=self.min_separation,
                                      exclude_border=self.exclude_border)
         else:
             xypos = self.xycoords

photutils/detection/irafstarfinder.py

@@ -40,8 +40,10 @@ class IRAFStarFinder(StarFinderBase):
         kernel in units of pixels.
 
     minsep_fwhm : float, optional
-        The minimum separation for detected objects in units of
-        ``fwhm``.
+        The separation (in units of ``fwhm``) for detected objects. The
+        minimum separation is calculated as ``int((fwhm * minsep_fwhm) +
+        0.5)`` and is clipped to a minimum value of 2. Note that large
+        values may result in long run times.
 
     sigma_radius : float, optional
         The truncation radius of the Gaussian kernel in units of sigma
@@ -87,11 +89,17 @@ class IRAFStarFinder(StarFinderBase):
             pixel values are negative. Therefore, setting ``peakmax`` to a
             non-positive value would result in exclusion of all objects.
 
-    xycoords : `None` or Nx2 `~numpy.ndarray`
+    xycoords : `None` or Nx2 `~numpy.ndarray`, optional
         The (x, y) pixel coordinates of the approximate centroid
         positions of identified sources. If ``xycoords`` are input, the
         algorithm will skip the source-finding step.
 
+    min_separation : `None` or float, optional
+        The minimum separation (in pixels) for detected objects. If
+        `None` then ``minsep_fwhm`` will be used, otherwise this keyword
+        overrides ``minsep_fwhm``. Note that large values may result in
+        long run times.
+
     Notes
     -----
     For the convolution step, this routine sets pixels beyond the image
@@ -131,7 +139,7 @@ class IRAFStarFinder(StarFinderBase):
     def __init__(self, threshold, fwhm, sigma_radius=1.5, minsep_fwhm=2.5,
                  sharplo=0.5, sharphi=2.0, roundlo=0.0, roundhi=0.2, sky=None,
                  exclude_border=False, brightest=None, peakmax=None,
-                 xycoords=None):
+                 xycoords=None, min_separation=None):
 
         if not np.isscalar(threshold):
             raise TypeError('threshold must be a scalar value.')
@@ -160,7 +168,14 @@ def __init__(self, threshold, fwhm, sigma_radius=1.5, minsep_fwhm=2.5,
 
         self.kernel = _StarFinderKernel(self.fwhm, ratio=1.0, theta=0.0,
                                         sigma_radius=self.sigma_radius)
-        self.min_separation = max(2, int((self.fwhm * self.minsep_fwhm) + 0.5))
+
+        if min_separation is not None:
+            if min_separation < 0:
+                raise ValueError('min_separation must be >= 0')
+            self.min_separation = min_separation
+        else:
+            self.min_separation = max(2, int((self.fwhm * self.minsep_fwhm)
+                                             + 0.5))
 
     @staticmethod
     def _validate_brightest(brightest):
@@ -472,8 +487,12 @@ def fwhm(self):
 
     @lazyproperty
     def roundness(self):
-        return np.sqrt(self.mu_diff**2
-                       + 4.0 * self.moments_central[:, 1, 1]**2) / self.mu_sum
+        # ignore divide-by-zero RuntimeWarning
+        with warnings.catch_warnings():
+            warnings.simplefilter('ignore', RuntimeWarning)
+            return (np.sqrt(self.mu_diff**2
+                            + 4.0 * self.moments_central[:, 1, 1]**2)
+                    / self.mu_sum)
 
     @lazyproperty
     def sharpness(self):

photutils/detection/starfinder.py

@@ -33,7 +33,8 @@ class StarFinder(StarFinderBase):
         A 2D array of the PSF kernel.
 
     min_separation : float, optional
-        The minimum separation for detected objects in pixels.
+        The minimum separation (in pixels) for detected objects. Note
+        that large values may result in long run times.
 
     exclude_border : bool, optional
         Whether to exclude sources found within half the size of the

photutils/detection/tests/test_daofinder.py

@@ -182,3 +182,17 @@ def test_invalid_xycoords(self):
         xycoords = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
         with pytest.raises(ValueError):
             DAOStarFinder(threshold=10, fwhm=1.5, xycoords=xycoords)
+
+    def test_min_separation(self):
+        threshold = 5
+        fwhm = 1.0
+        finder1 = DAOStarFinder(threshold, fwhm, sigma_radius=1.5)
+        tbl1 = finder1(DATA)
+        finder2 = DAOStarFinder(threshold, fwhm, sigma_radius=1.5,
+                                min_separation=3.0)
+        tbl2 = finder2(DATA)
+        assert len(tbl1) > len(tbl2)
+
+        match = 'min_separation must be >= 0'
+        with pytest.raises(ValueError, match=match):
+            DAOStarFinder(threshold=10, fwhm=1.5, min_separation=-1.0)

photutils/detection/tests/test_irafstarfinder.py

@@ -150,3 +150,22 @@ def test_invalid_xycoords(self):
         xycoords = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
         with pytest.raises(ValueError):
             IRAFStarFinder(threshold=10, fwhm=1.5, xycoords=xycoords)
+
+    def test_min_separation(self):
+        threshold = 5
+        fwhm = 1.0
+        starfinder1 = IRAFStarFinder(threshold, fwhm, sigma_radius=1.5)
+        tbl1 = starfinder1(DATA)
+        starfinder2 = IRAFStarFinder(threshold, fwhm, sigma_radius=1.5,
+                                     min_separation=3.0)
+        tbl2 = starfinder2(DATA)
+        assert np.all(tbl1 == tbl2)
+
+        starfinder3 = IRAFStarFinder(threshold, fwhm, sigma_radius=1.5,
+                                     min_separation=2.0)
+        tbl3 = starfinder3(DATA)
+        assert len(tbl3) > len(tbl2)
+
+        match = 'min_separation must be >= 0'
+        with pytest.raises(ValueError, match=match):
+            IRAFStarFinder(threshold=10, fwhm=1.5, min_separation=-1.0)