Repository: astropy/photutils

Issue:
Rename 'flux' column in tables generated by ``StarFinderBase`` objects
I think it makes sense to rename the ``flux`` column in tables generated by ``StarFinderBase`` objects to something like ``flux_peak_density``. I'm saying this because the name ``flux`` induces users to use those numbers as the initial guesses for the flux parameter in the PSF photometry code.

P.S.: I myself used ``flux`` from those tables as initial guesses for PSF photometry which end up generating bad results because those values were too far from the "true" flux value, thanks @eteq saw that and we end up using aperture photometry to generate quite good guesses.

Discussion:

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 4a94f8b7..31501f9c 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -163,6 +163,14 @@ API Changes
     ``IRAFStarFinder`` has been removed. Also, there will no longer be a
     ``sky`` column in the output table. [#1884]
 
+  - The ``DAOStarFinder`` ``flux`` and ``mag`` columns were changed
+    to give sensible values. Previously, the ``flux`` value was defined
+    by the original DAOFIND algorithm as a measure of the intensity
+    ratio of the amplitude of the best fitting Gaussian function at the
+    object position to the detection threshold. A ``daofind_mag`` column
+    was added for comparison to the original IRAF DAOFIND algorithm.
+    [#1885]
+
 - ``photutils.isophote``
 
   - The ``build_ellipse_model`` function now raises a ``ValueError`` if
diff --git a/docs/detection.rst b/docs/detection.rst
index 7da0bb04..d1bcba39 100644
--- a/docs/detection.rst
+++ b/docs/detection.rst
@@ -75,19 +75,19 @@ finder::
     ...     if col not in ('id', 'npix'):
     ...         sources[col].info.format = '%.2f'  # for consistent table output
     >>> sources.pprint(max_width=76)  # doctest: +REMOTE_DATA, +FLOAT_CMP
-     id xcentroid ycentroid sharpness roundness1 ... npix   peak  flux  mag
-    --- --------- --------- --------- ---------- ... ---- ------- ---- -----
-      1    144.25      6.38      0.58       0.20 ...   25 6903.00 5.70 -1.89
-      2    208.67      6.82      0.48      -0.13 ...   25 7896.00 6.72 -2.07
-      3    216.93      6.58      0.69      -0.71 ...   25 2195.00 1.67 -0.55
-      4    351.63      8.55      0.49      -0.34 ...   25 6977.00 5.90 -1.93
-      5    377.52     12.07      0.52       0.37 ...   25 1260.00 1.12 -0.12
-    ...       ...       ...       ...        ... ...  ...     ...  ...   ...
-    282    267.90    398.62      0.27      -0.43 ...   25 9299.00 5.44 -1.84
-    283    271.47    398.91      0.37       0.19 ...   25 8028.00 5.07 -1.76
-    284    299.05    398.78      0.26      -0.67 ...   25 9072.00 5.56 -1.86
-    285    299.99    398.77      0.29       0.36 ...   25 9253.00 5.32 -1.82
-    286    360.45    399.52      0.37      -0.19 ...   25 8079.00 6.92 -2.10
+     id xcentroid ycentroid sharpness ...   peak     flux    mag   daofind_mag
+    --- --------- --------- --------- ... ------- --------- ------ -----------
+      1    144.25      6.38      0.58 ... 6903.00  45735.00 -11.65       -1.89
+      2    208.67      6.82      0.48 ... 7896.00  62118.00 -11.98       -2.07
+      3    216.93      6.58      0.69 ... 2195.00  12436.00 -10.24       -0.55
+      4    351.63      8.55      0.49 ... 6977.00  55313.00 -11.86       -1.93
+      5    377.52     12.07      0.52 ... 1260.00   9078.00  -9.89       -0.12
+    ...       ...       ...       ... ...     ...       ...    ...         ...
+    282    267.90    398.62      0.27 ... 9299.00 147372.00 -12.92       -1.84
+    283    271.47    398.91      0.37 ... 8028.00 115913.00 -12.66       -1.76
+    284    299.05    398.78      0.26 ... 9072.00 140781.00 -12.87       -1.86
+    285    299.99    398.77      0.29 ... 9253.00 142233.00 -12.88       -1.82
+    286    360.45    399.52      0.37 ... 8079.00  81455.00 -12.28       -2.10
     Length = 286 rows
 
 Let's plot the image and mark the location of detected sources:
diff --git a/docs/getting_started.rst b/docs/getting_started.rst
index 3fb2a0b2..56fdcb80 100644
--- a/docs/getting_started.rst
+++ b/docs/getting_started.rst
@@ -34,21 +34,21 @@ detected sources are returned as an Astropy `~astropy.table.Table`::
     >>> daofind = DAOStarFinder(fwhm=4.0, threshold=3.0 * bkg_sigma)  # doctest: +REMOTE_DATA
     >>> sources = daofind(image)  # doctest: +REMOTE_DATA
     >>> for col in sources.colnames:  # doctest: +REMOTE_DATA
-    ...     sources[col].info.format = '%.8g'  # for consistent table output
+    ...     sources[col].info.format = '%.5g'  # for consistent table output
     >>> print(sources)  # doctest: +REMOTE_DATA
-     id xcentroid ycentroid  sharpness  ... npix peak    flux       mag
-    --- --------- ---------- ---------- ... ---- ---- --------- -----------
-      1 182.83866 0.16767019 0.85099873 ...   25 3824 2.8028346  -1.1189937
-      2 189.20431 0.26081353  0.7400477 ...   25 4913 3.8729185  -1.4700959
-      3 5.7946491  2.6125424 0.39589731 ...   25 7752 4.1029107  -1.5327302
-      4 36.847063  1.3220228 0.29594528 ...   25 8739 7.4315818  -2.1777032
-      5 3.2565602   5.418952 0.35985495 ...   25 6935 3.8126298  -1.4530616
-    ...       ...        ...        ... ...  ...  ...       ...         ...
-    148 124.31327  188.30523  0.5362742 ...   25 6702 6.6358543  -2.0547421
-    149 24.257207  194.71494 0.44169546 ...   25 8342 3.2671037  -1.2854073
-    150    116.45  195.05923 0.67080547 ...   25 3299 2.8775221  -1.1475467
-    151 18.958086  196.34207 0.56502139 ...   25 3854 2.3835296 -0.94305138
-    152 111.52575  195.73192 0.45827852 ...   25 8109 7.9278607    -2.24789
+     id xcentroid ycentroid sharpness ... peak    flux      mag   daofind_mag
+    --- --------- --------- --------- ... ---- ---------- ------- -----------
+      1    182.84   0.16767     0.851 ... 3824      15750 -10.493      -1.119
+      2     189.2   0.26081   0.74005 ... 4913      25870 -11.032     -1.4701
+      3    5.7946    2.6125    0.3959 ... 7752 1.3525e+05 -12.828     -1.5327
+      4    36.847     1.322   0.29595 ... 8739 1.3158e+05 -12.798     -2.1777
+      5    3.2566     5.419   0.35985 ... 6935 1.2454e+05 -12.738     -1.4531
+    ...       ...       ...       ... ...  ...        ...     ...         ...
+    148    124.31    188.31   0.53627 ... 6702      44592 -11.623     -2.0547
+    149    24.257    194.71    0.4417 ... 8342 1.5645e+05 -12.986     -1.2854
+    150    116.45    195.06   0.67081 ... 3299      17418 -10.602     -1.1475
+    151    18.958    196.34   0.56502 ... 3854      48800 -11.721    -0.94305
+    152    111.53    195.73   0.45828 ... 8109      74247 -12.177     -2.2479
     Length = 152 rows
 
 Using the source locations (i.e., the ``xcentroid`` and ``ycentroid``
diff --git a/docs/psf.rst b/docs/psf.rst
index f6feee5a..94e44091 100644
--- a/docs/psf.rst
+++ b/docs/psf.rst
@@ -478,19 +478,20 @@ astropy table)::
     >>> psfphot.finder_results['peak'].info.format = '.4f'
     >>> psfphot.finder_results['flux'].info.format = '.4f'
     >>> psfphot.finder_results['mag'].info.format = '.4f'
+    >>> psfphot.finder_results['daofind_mag'].info.format = '.4f'
     >>> print(psfphot.finder_results)  # doctest: +FLOAT_CMP
-     id xcentroid ycentroid sharpness ... npix   peak   flux    mag
-    --- --------- --------- --------- ... ---- ------- ------ -------
-      1   54.5299    7.7460    0.6006 ...   25 53.5953 6.9777 -2.1093
-      2   29.0927   25.5992    0.5955 ...   25 57.1982 7.5568 -2.1958
-      3   79.6185   28.7515    0.5957 ...   25 65.7175 8.6306 -2.3401
-      4   63.2485   48.6134    0.5802 ...   25 58.3985 7.7333 -2.2209
-      5   88.8820   54.1311    0.5948 ...   25 69.1869 9.4443 -2.4379
-      6   79.8727   61.1208    0.6216 ...   25 74.0935 9.8165 -2.4799
-      7   90.9621   72.0803    0.6167 ...   25 68.4157 9.1492 -2.4035
-      8    7.7962   78.5465    0.5979 ...   25 66.2173 8.4468 -2.3167
-      9    5.5858   89.8664    0.5741 ...   25 54.3786 7.0389 -2.1188
-     10   71.8303   90.5624    0.6038 ...   25 73.5747 9.5639 -2.4516
+     id xcentroid ycentroid sharpness ...   peak    flux     mag   daofind_mag
+    --- --------- --------- --------- ... ------- -------- ------- -----------
+      1   54.5299    7.7460    0.6006 ... 53.5953 476.3221 -6.6948     -2.1093
+      2   29.0927   25.5992    0.5955 ... 57.1982 499.4443 -6.7462     -2.1958
+      3   79.6185   28.7515    0.5957 ... 65.7175 574.1382 -6.8975     -2.3401
+      4   63.2485   48.6134    0.5802 ... 58.3985 521.4656 -6.7931     -2.2209
+      5   88.8820   54.1311    0.5948 ... 69.1869 576.2842 -6.9016     -2.4379
+      6   79.8727   61.1208    0.6216 ... 74.0935 612.8353 -6.9684     -2.4799
+      7   90.9621   72.0803    0.6167 ... 68.4157 561.7163 -6.8738     -2.4035
+      8    7.7962   78.5465    0.5979 ... 66.2173 595.6881 -6.9375     -2.3167
+      9    5.5858   89.8664    0.5741 ... 54.3786 505.6093 -6.7595     -2.1188
+     10   71.8303   90.5624    0.6038 ... 73.5747 639.9299 -7.0153     -2.4516
 
 The ``fit_info`` attribute contains a dictionary with detailed
 information returned from the ``fitter`` for each source::
diff --git a/photutils/detection/core.py b/photutils/detection/core.py
index 2961f36b..446d70c9 100644
--- a/photutils/detection/core.py
+++ b/photutils/detection/core.py
@@ -237,6 +237,7 @@ class _StarFinderKernel:
         denom = (self.a * self.c) - self.b**2
 
         # nx and ny are always odd
+        # minimum kernel size is 5x5
         self.nx = 2 * int(max(2, math.sqrt(self.c * self.f / denom))) + 1
         self.ny = 2 * int(max(2, math.sqrt(self.a * self.f / denom))) + 1
 
diff --git a/photutils/detection/daofinder.py b/photutils/detection/daofinder.py
index 61266d8b..1bdb2a3a 100644
--- a/photutils/detection/daofinder.py
+++ b/photutils/detection/daofinder.py
@@ -247,7 +247,8 @@ class DAOStarFinder(StarFinderBase):
         Returns
         -------
         table : `~astropy.table.QTable` or `None`
-            A table of found stars with the following parameters:
+            A table of found stars. `None` is returned if no stars are
+            found. The table contains the following parameters:
 
             * ``id``: unique object identification number.
             * ``xcentroid, ycentroid``: object centroid.
@@ -258,26 +259,17 @@ class DAOStarFinder(StarFinderBase):
             * ``npix``: the total number of pixels in the Gaussian kernel
               array.
             * ``peak``: the peak pixel value of the object.
-            * ``flux``: the object DAOFIND "flux" calculated as the peak
-              density in the convolved image divided by the detection
-              threshold. This derivation matches that of DAOFIND.
+            * ``flux``: the object instrumental flux calculated as the
+              sum of data values within the kernel footprint.
             * ``mag``: the object instrumental magnitude calculated as
-              ``-2.5 * log10(flux)``. The derivation matches that of
-              DAOFIND.
-
-            `None` is returned if no stars are found.
-
-            .. warning::
-                The ``flux`` parameter returned by the DAOFIND algorithm
-                is a bit of a misnomer. It is not a typical integrated
-                source flux, and its value critically depends on the
-                input ``threshold`` parameter. The ``flux`` parameter is
-                calculated as the peak density in the convolved image
-                divided by the effective detection threshold. ``flux``
-                and the corresponding ``mag`` parameter are reported
-                in the output table for comparison to the IRAF DAOFIND
-                output only. They are not intended to be used as a
-                measure of source flux.
+              ``-2.5 * log10(flux)``.
+            * ``daofind_mag``: the "mag" parameter returned by the DAOFIND
+              algorithm. It is a measure of the intensity ratio of the
+              amplitude of the best fitting Gaussian function at the
+              object position to the detection threshold. This parameter
+              is reported only for comparison to the IRAF DAOFIND
+              output. It should not be interpreted as a magnitude
+              derived from an integrated flux.
         """
         # here we validate the units, but do not strip them
         inputs = (data, self.threshold, self.peakmax)
@@ -382,7 +374,7 @@ class _DAOStarFinderCatalog:
         self.cutout_center = tuple((size - 1) // 2 for size in kernel.shape)
         self.default_columns = ('id', 'xcentroid', 'ycentroid', 'sharpness',
                                 'roundness1', 'roundness2', 'npix', 'peak',
-                                'flux', 'mag')
+                                'flux', 'mag', 'daofind_mag')
 
     def __len__(self):
         return len(self.xypos)
@@ -691,22 +683,12 @@ class _DAOStarFinderCatalog:
 
     @lazyproperty
     def flux(self):
-        """
-        The peak density in the convolved image divided by the effective
-        detection threshold.
-
-        This is DAOStarFinder's definition of "flux" and matches that of
-        DAOFIND.
-        """
-        # ignore divide-by-zero RuntimeWarning if threshold = 0
-        with warnings.catch_warnings():
-            warnings.simplefilter('ignore', RuntimeWarning)
-            flux = self.convdata_peak / self.threshold_eff
-
+        fluxes = [np.sum(arr) for arr in self.cutout_data]
         if self.unit is not None:
-            flux = flux.value * self.unit
-
-        return flux
+            fluxes = u.Quantity(fluxes)
+        else:
+            fluxes = np.array(fluxes)
+        return fluxes
 
     @lazyproperty
     def mag(self):
@@ -718,6 +700,20 @@ class _DAOStarFinderCatalog:
                 flux = flux.value
             return -2.5 * np.log10(flux)
 
+    @lazyproperty
+    def daofind_mag(self):
+        """
+        The "mag" parameter returned by the original DAOFIND algorithm.
+
+        It is a measure of the intensity ratio of the amplitude of the
+        best fitting Gaussian function at the object position to the
+        detection threshold.
+        """
+        # ignore RunTimeWarning if flux is <= 0
+        with warnings.catch_warnings():
+            warnings.simplefilter('ignore', category=RuntimeWarning)
+            return -2.5 * np.log10(self.convdata_peak / self.threshold_eff)
+
     @lazyproperty
     def npix(self):
         return np.full(len(self), fill_value=self.kernel.data.size)
diff --git a/photutils/detection/irafstarfinder.py b/photutils/detection/irafstarfinder.py
index 779c346c..bb2ddfb2 100644
--- a/photutils/detection/irafstarfinder.py
+++ b/photutils/detection/irafstarfinder.py
@@ -116,20 +116,23 @@ class IRAFStarFinder(StarFinderBase):
     input parameters. The equivalent input values for `IRAFStarFinder`
     are:
 
-    * ``fwhm = hwhmpsf * 2``
-    * ``sigma_radius = fradius * sqrt(2.0*log(2.0))``
-    * ``minsep_fwhm = 0.5 * sepmin``
+      * ``fwhm = hwhmpsf * 2``
+      * ``sigma_radius = fradius * sqrt(2.0*log(2.0))``
+      * ``minsep_fwhm = 0.5 * sepmin``
 
     The main differences between `~photutils.detection.DAOStarFinder`
     and `~photutils.detection.IRAFStarFinder` are:
 
-    * `~photutils.detection.IRAFStarFinder` always uses a 2D
-      circular Gaussian kernel, while
-      `~photutils.detection.DAOStarFinder` can use an elliptical
-      Gaussian kernel.
+      * `~photutils.detection.IRAFStarFinder` always uses a 2D
+        circular Gaussian kernel, while
+        `~photutils.detection.DAOStarFinder` can use an elliptical
+        Gaussian kernel.
 
-    * `~photutils.detection.IRAFStarFinder` calculates the objects'
-      centroid, roundness, and sharpness using image moments.
+      * `IRAFStarFinder` internally calculates a "sky" background level
+        based on unmasked pixels within the kernel footprint.
+
+      * `~photutils.detection.IRAFStarFinder` calculates the objects'
+        centroid, roundness, and sharpness using image moments.
     """
 
     def __init__(self, threshold, fwhm, sigma_radius=1.5, minsep_fwhm=2.5,
@@ -222,7 +225,8 @@ class IRAFStarFinder(StarFinderBase):
         Returns
         -------
         table : `~astropy.table.QTable` or `None`
-            A table of found objects with the following parameters:
+            A table of found stars. `None` is returned if no stars are
+            found. The table contains the following parameters:
 
             * ``id``: unique object identification number.
             * ``xcentroid, ycentroid``: object centroid.
@@ -233,11 +237,11 @@ class IRAFStarFinder(StarFinderBase):
               the positive x axis).
             * ``npix``: the total number of (positive) unmasked pixels.
             * ``peak``: the peak, sky-subtracted, pixel value of the object.
-            * ``flux``: the object instrumental flux.
+            * ``flux``: the object instrumental flux calculated as the
+               sum of sky-subtracted data values within the kernel
+               footprint.
             * ``mag``: the object instrumental magnitude calculated as
               ``-2.5 * log10(flux)``.
-
-            `None` is returned if no stars are found.
         """
         inputs = (data, self.threshold, self.peakmax)
         names = ('data', 'threshold', 'peakmax')


Test patch:
diff --git a/photutils/detection/tests/test_daofinder.py b/photutils/detection/tests/test_daofinder.py
index 70f518e8..570d6cb9 100644
--- a/photutils/detection/tests/test_daofinder.py
+++ b/photutils/detection/tests/test_daofinder.py
@@ -24,12 +24,12 @@ class TestDAOStarFinder:
         tbl1 = finder1(data << units)
         assert_array_equal(tbl0, tbl1)
 
+        assert np.min(tbl0['flux']) > 150
+
         # test that sources are returned with threshold = 0
         finder = DAOStarFinder(0, fwhm)
         tbl = finder(data)
         assert len(tbl) == 25
-        assert np.all(~np.isfinite(tbl['flux']))
-        assert np.all(~np.isfinite(tbl['mag']))
 
     def test_daofind_inputs(self):
         match = 'threshold must be a scalar value'