Merge pull request #240 from larrybradley/find-peaks-table

Return an Astropy Table from find_peaks()

CHANGES.rst

@@ -1,4 +1,14 @@
 0.2 (unreleased)
 ----------------
 
-- No changes yet
+General
+^^^^^^^
+
+
+New Features
+^^^^^^^^^^^^
+
+- ``photutils.detection``
+
+  - ``find_peaks`` now returns an Astropy Table containing the (x, y)
+    positions and peak values. [#240]

docs/photutils/detection.rst

@@ -253,7 +253,20 @@ sigma above the background and a separated by a least 2 pixels:
     >>> data = make_100gaussians_image()
     >>> mean, median, std = sigma_clipped_stats(data, sigma=3.0)
     >>> threshold = median + (10.0 * std)
-    >>> peaks = find_peaks(data, threshold, min_separation=2)
+    >>> tbl = find_peaks(data, threshold, min_separation=2)
+    >>> print(tbl[:10])    # print only the first 10 peaks
+    x_peak y_peak   peak_value
+    ------ ------ -------------
+       289     22 35.8532759965
+       442     31 30.2399941373
+        89     59 41.2190469279
+         7     70 33.2880647048
+       258     75 26.5624808518
+       463     80 28.7588206692
+       182     93 38.0885687202
+       227    104  29.291381982
+        71    112 35.0913969385
+       206    114 29.7546446436
 
 And let's plot the location of the detected peaks in the image:
 
@@ -264,8 +277,8 @@ And let's plot the location of the detected peaks in the image:
     >>> import matplotlib.pylab as plt
     >>> norm = ImageNormalize(stretch=SqrtStretch())
     >>> plt.imshow(data, cmap='Greys_r', origin='lower', norm=norm)
-    >>> plt.plot(peaks.T[1], peaks.T[0], ls='none', color='cyan', marker='+',
-    ...          ms=10, lw=1.5)
+    >>> plt.plot(tbl['x_peak'], tbl['y_peak'], ls='none', color='cyan',
+    ...          marker='+', ms=10, lw=1.5)
     >>> plt.xlim(0, data.shape[1]-1)
     >>> plt.ylim(0, data.shape[0]-1)
 
@@ -277,15 +290,15 @@ And let's plot the location of the detected peaks in the image:
     data = make_100gaussians_image()
     mean, median, std = sigma_clipped_stats(data, sigma=3.0)
     threshold = median + (10.0 * std)
-    peaks = find_peaks(data, threshold, min_separation=2)
+    tbl = find_peaks(data, threshold, min_separation=2)
 
     from photutils.extern.imageutils.visualization import SqrtStretch
     from photutils.extern.imageutils.visualization.mpl_normalize import ImageNormalize
     import matplotlib.pylab as plt
     norm = ImageNormalize(stretch=SqrtStretch())
     plt.imshow(data, cmap='Greys_r', origin='lower', norm=norm)
-    plt.plot(peaks.T[1], peaks.T[0], ls='none', color='cyan', marker='+',
-             ms=10, lw=1.5)
+    plt.plot(tbl['x_peak'], tbl['y_peak'], ls='none', color='cyan',
+             marker='+', ms=10, lw=1.5)
     plt.xlim(0, data.shape[1]-1)
     plt.ylim(0, data.shape[0]-1)
 

photutils/detection/core.py

@@ -4,6 +4,7 @@
 from __future__ import (absolute_import, division, print_function,
                         unicode_literals)
 import numpy as np
+from astropy.table import Table
 from astropy.convolution import Kernel2D
 from ..extern.imageutils.stats import sigma_clipped_stats
 
@@ -301,9 +302,9 @@ def find_peaks(data, threshold, min_separation=2, exclude_border=True,
 
     Returns
     -------
-    output : `~numpy.ndarray`
-        An ``Nx2`` array where the rows contain the ``(y, x)`` pixel
-        coordinates of the local peaks.
+    output : `~astropy.table.Table`
+        A table containing the x and y pixel location of the peaks and
+        their values.
     """
 
     if segment_image is not None:
@@ -317,10 +318,15 @@ def find_peaks(data, threshold, min_separation=2, exclude_border=True,
                             exclude_border=exclude_border, indices=True,
                             num_peaks=npeaks, footprint=footprint,
                             labels=segment_image)
-    if coords.shape[0] <= npeaks:
-        return coords
-    else:
+    y_peaks, x_peaks = coords[:, 0], coords[:, 1]
+    peak_values = data[y_peaks, x_peaks]
+
+    if coords.shape[0] > npeaks:
         # NOTE: num_peaks is ignored by peak_local_max() if labels are input
-        peak_values = data[coords[:, 0], coords[:, 1]]
-        idx_maxsort = np.argsort(peak_values)[::-1]
-        return coords[idx_maxsort][:npeaks]
+        idx = np.argsort(peak_values)[::-1][:npeaks]
+        x_peaks = x_peaks[idx]
+        y_peaks = y_peaks[idx]
+        peak_values = peak_values[idx]
+
+    return Table((x_peaks, y_peaks, peak_values),
+                 names=('x_peak', 'y_peak', 'peak_value'))

photutils/detection/findstars.py

@@ -16,7 +16,6 @@
 __all__ = ['daofind', 'irafstarfind']
 
 
-
 def daofind(data, 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):
@@ -389,9 +388,9 @@ def _findobjs(data, threshold, kernel, min_separation=None,
         else:
             from skimage.morphology import disk
             footprint = disk(min_separation)
-        coords = find_peaks(convolved_data, threshold,
-                            exclude_border=exclude_border,
-                            footprint=footprint)
+        tbl = find_peaks(convolved_data, threshold,
+                         exclude_border=exclude_border, footprint=footprint)
+        coords = np.transpose([tbl['y_peak'], tbl['x_peak']])
     else:
         object_slices = ndimage.find_objects(object_labels)
         coords = []

photutils/detection/tests/test_core.py

@@ -191,22 +191,25 @@ def test_filtering_kernel(self):
 class TestFindPeaks(object):
     def test_find_peaks(self):
         """Test basic peak detection."""
-        coords = find_peaks(PEAKDATA, 0.1, min_separation=1,
-                            exclude_border=False)
-        assert_array_equal(coords, PEAKREF1)
+        tbl = find_peaks(PEAKDATA, 0.1, min_separation=1,
+                         exclude_border=False)
+        assert_array_equal(tbl['x_peak'], PEAKREF1[:, 1])
+        assert_array_equal(tbl['y_peak'], PEAKREF1[:, 0])
+        assert_array_equal(tbl['peak_value'], [1., 1.])
 
     def test_segment_image(self):
         segm = PEAKDATA.copy()
-        coords = find_peaks(PEAKDATA, 0.1, min_separation=1,
-                            exclude_border=False, segment_image=segm)
-        assert_array_equal(coords, PEAKREF1)
+        tbl = find_peaks(PEAKDATA, 0.1, min_separation=1,
+                         exclude_border=False, segment_image=segm)
+        assert_array_equal(tbl['x_peak'], PEAKREF1[:, 1])
+        assert_array_equal(tbl['y_peak'], PEAKREF1[:, 0])
 
     def test_segment_image_npeaks(self):
         segm = PEAKDATA.copy()
-        coords = find_peaks(PEAKDATA, 0.1, min_separation=1,
-                            exclude_border=False, segment_image=segm,
-                            npeaks=1)
-        assert_array_equal(coords, np.array([PEAKREF1[1]]))
+        tbl = find_peaks(PEAKDATA, 0.1, min_separation=1,
+                         exclude_border=False, segment_image=segm, npeaks=1)
+        assert_array_equal(tbl['x_peak'], PEAKREF1[1, 1])
+        assert_array_equal(tbl['y_peak'], PEAKREF1[1, 0])
 
     def test_segment_image_shape(self):
         segm = np.zeros((2, 2))
@@ -215,17 +218,15 @@ def test_segment_image_shape(self):
 
     def test_exclude_border(self):
         """Test exclude_border."""
-        coords = find_peaks(PEAKDATA, 0.1, min_separation=1,
-                            exclude_border=True)
-        assert_array_equal(coords, PEAKREF2)
+        tbl = find_peaks(PEAKDATA, 0.1, min_separation=1, exclude_border=True)
+        assert_array_equal(len(tbl), 0)
 
     def test_zerodet(self):
         """Test with large threshold giving no sources."""
-        coords = find_peaks(PEAKDATA, 5., min_separation=1,
-                            exclude_border=True)
-        assert_array_equal(coords, PEAKREF2)
+        tbl = find_peaks(PEAKDATA, 5., min_separation=1, exclude_border=True)
+        assert_array_equal(len(tbl), 0)
 
     def test_min_separation_int(self):
-        coords1 = find_peaks(PEAKDATA, 0.1, min_separation=2)
-        coords2 = find_peaks(PEAKDATA, 0.1, min_separation=2.5)
-        assert_array_equal(coords1, coords2)
+        tbl1 = find_peaks(PEAKDATA, 0.1, min_separation=2.)
+        tbl2 = find_peaks(PEAKDATA, 0.1, min_separation=2.5)
+        assert_array_equal(tbl1, tbl2)