Merge pull request #1666 from larrybradley/depth-coords

Improve random apertures in ImageDepth

CHANGES.rst

@@ -64,6 +64,11 @@ New Features
   - Added a ``reset_cmap`` method to ``SegmentationImage`` for resetting
     the colormap to a new random colormap. [#1649]
 
+- ``photutils.utils``
+
+  - Improved the generation of random aperture positions in
+    ``ImageDepth``. [#1666]
+
 Bug Fixes
 ^^^^^^^^^
 

photutils/utils/depths.py

@@ -10,6 +10,7 @@
 from astropy.stats import SigmaClip
 from astropy.utils.exceptions import AstropyUserWarning
 
+from photutils.utils._coords import apply_separation
 from photutils.utils._progress_bars import add_progress_bar
 from photutils.utils.footprints import circular_footprint
 
@@ -143,7 +144,7 @@ class ImageDepth:
     ...                    progress_bar=False)
     >>> limits = depth(data, mask)
     >>> print(limits)  # doctest: +FLOAT_CMP
-    (67.24989263150033, 19.330771011755676)
+    (68.0112578151062, 19.318547982855563)
 
     .. plot::
         :include-source:
@@ -270,7 +271,11 @@ def __call__(self, data, mask):
             if self.overlap:
                 xycoords = self._make_coords(all_xycoords, napers)
             else:
-                xycoords = self._make_nonoverlap_coords(all_xycoords, napers)
+                # cut the number of coords (only need to input ~10x)
+                xycoords = self._make_coords(all_xycoords, napers * 10)
+                min_separation = self.aper_radius * 2.0
+                xycoords = apply_separation(xycoords, min_separation)
+                xycoords = xycoords[0:self.napers]
 
             apers = CircularAperture(xycoords, r=self.aper_radius)
             apertures.append(apers)
@@ -509,67 +514,3 @@ def _make_coords(self, xycoords, napers):
         xycoords += shift
 
         return xycoords
-
-    def _make_nonoverlap_coords(self, init_xycoords, napers):
-        """
-        Randomly choose ``napers`` (without replacement) coordinates
-        from the input ``xycoords`` that do not overlap.
-
-        Parameters
-        ----------
-        xycoords : 2xN `~numpy.ndarray`
-            The (x, y) coordinates.
-        napers : int
-            The number of aperture to make.
-
-        Returns
-        --------
-        xycoords : 2xN `~numpy.ndarray`
-            The (x, y) coordinates.
-        """
-        from scipy.spatial import KDTree
-
-        minsep = self.aper_radius * 2.0
-        xycoords = np.zeros((0, 2))  # placeholder for while loop
-
-        # attempt to generate all the coordinates at once; this will
-        # work only for non-crowded blank areas
-        niter = 1
-        while xycoords.shape[0] < self.napers:
-            if niter > 10:
-                break
-
-            new_xycoords = self._make_coords(init_xycoords, napers)
-            if niter == 1:
-                xycoords = new_xycoords
-            else:
-                xycoords = np.vstack((xycoords, new_xycoords))
-
-            dist, _ = KDTree(xycoords).query(xycoords, k=[2])
-            mask = (dist >= minsep).squeeze()
-            xycoords = xycoords[mask]
-            niter += 1
-
-        # add new coordinates one by one (slower)
-        napers = len(init_xycoords)
-        if xycoords.shape[0] < self.napers:
-            new_xycoords = self._make_coords(init_xycoords, napers)
-            if xycoords.shape[0] == 0:
-                xycoords = new_xycoords[0]
-                new_xycoords = new_xycoords[1:]
-
-            count = 1
-            while (xycoords.shape[0] < self.napers
-                   and count < self.overlap_maxiters):
-                idx = self.rng.choice(new_xycoords.shape[0], 1, replace=False)
-                new_xy = new_xycoords[idx, :]
-                dist, _ = KDTree(new_xy).query(xycoords, 1)
-                if np.min(dist) > minsep:
-                    xycoords = np.vstack((xycoords, new_xy))
-                    count = 0
-                else:
-                    count += 1
-
-        xycoords = xycoords[0:self.napers]
-
-        return xycoords

photutils/utils/tests/test_depths.py

@@ -48,7 +48,7 @@ def test_image_depth(self, units, overlap):
         if overlap:
             exp_limits = (66.60324687694235, 19.341261496655026)
         else:
-            exp_limits = (67.65345151009167, 19.324275104703975)
+            exp_limits = (68.60472024061636, 19.30911500577945)
 
         data = self.data
         fluxlim = exp_limits[0]
@@ -82,7 +82,7 @@ def test_many_apertures(self):
         with pytest.raises(ValueError):
             depth(self.data, mask)
 
-        depth = ImageDepth(radius, nsigma=5.0, napers=200, niters=2,
+        depth = ImageDepth(radius, nsigma=5.0, napers=250, niters=2,
                            overlap=False, seed=123, zeropoint=23.9,
                            progress_bar=False)
         mask = np.zeros(self.data.shape)