No more centroiding in Simple Aperture Photometry plugin (#1841)

* No more centroiding in
Simple Aperture Photometry plugin.

Removed outdated known issue.

Updated docs and tests.

* DOC: Clarify how to recenter

* DOC: Fix grammar

Co-authored-by: Ricky O'Steen <39831871+rosteen@users.noreply.github.com>

Author: pllim

CHANGES.rst

@@ -54,6 +54,11 @@ Cubeviz
 Imviz
 ^^^^^
 
+- Simple Aperture Photometry plugin no longer performs centroiding.
+  For radial profile, curve of growth, and table reporting, the aperture
+  center is used instead. For centroiding, use "Recenter" feature in
+  the Subset Tools plugin. [#1841]
+
 Mosviz
 ^^^^^^
 

docs/imviz/export_data.rst

@@ -51,11 +51,8 @@ The columns are as follow:
 
 * :attr:`~photutils.aperture.ApertureStats.id`: ID number assigned to the row,
   starting from 1.
-* :attr:`~photutils.aperture.ApertureStats.xcentroid`,
-  :attr:`~photutils.aperture.ApertureStats.ycentroid`: Pixel centroids
-  calculated using moments. This might differ from center of the aperture.
-* :attr:`~photutils.aperture.ApertureStats.sky_centroid`:
-  `~astropy.coordinates.SkyCoord` associated with the centroid.
+* ``xcenter``, ``ycenter``: Center of the aperture (0-indexed).
+* ``sky_center``: `~astropy.coordinates.SkyCoord` associated with the center.
   If WCS is not available, this field is `None`.
 * ``background``: The value from :guilabel:`Background value`, with unit attached.
 * :attr:`~photutils.aperture.ApertureStats.sum`: Sum of flux in the aperture.
@@ -108,9 +105,8 @@ The columns are as follow:
 .. note::
 
     Aperture sum and statistics are done on the originally drawn aperture only.
-    Even though centroid is calculated, it is not used to move the aperture
-    to the new center. However, radial profiles (including Gaussian fitting, if any)
-    and curve of growth do use the centroid as zero-point on the X-axis.
+    You can use the :ref:`imviz-subset-plugin` plugin to center it first on the
+    object of interest, if you wish.
 
 Once you have the results in a table, you can further manipulated them as
 documented in :ref:`astropy:astropy-table`.

docs/imviz/plugins.rst

@@ -143,6 +143,10 @@ an interactively selected region. A typical workflow is as follows:
 2. Draw a region over the object of interest (see :ref:`imviz_defining_spatial_regions`).
 3. Select the desired image using the :guilabel:`Data` dropdown menu.
 4. Select the desired region using the :guilabel:`Subset` dropdown menu.
+   You can use the :ref:`imviz-subset-plugin` plugin to center it first on the
+   object of interest using its center of mass, if you wish.
+   Depending on the object, it may take several iterations for re-centering
+   to converge, or it may never converge at all.
 5. If you want to subtract background before performing photometry,
    you have the following 3 options. Otherwise if your image is already
    background subtracted, choose "Manual" and leave the background set at 0:
@@ -240,7 +244,7 @@ catalog dropdown menu.
     This plugin is still under active development. As a result, the search only uses the SDSS DR17 catalog
     and works best when you only have a single image loaded in a viewer.
 
-To load a catalog from a supported `JWST ECSV catalog file <https://jwst-pipeline.readthedocs.io/en/latest/jwst/source_catalog/main.html#output-products>`_, choose "From File...".  
+To load a catalog from a supported `JWST ECSV catalog file <https://jwst-pipeline.readthedocs.io/en/latest/jwst/source_catalog/main.html#output-products>`_, choose "From File...".
 The file must be able to be parsed by `astropy.table.Table.read` and contain a column labeled 'sky_centroid'.
 Clicking :guilabel:`SEARCH` will show markers for any entry within the filtered zoom window.
 

docs/known_bugs.rst

@@ -143,14 +143,6 @@ to show the markers, or not at all. This is a known bug reported in
 https://github.com/glue-viz/glue-jupyter/issues/243 . If you encounter this,
 try a different OS/browser combo.
 
-Simple Aperture Photometry Plugin and dithered images
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-Due to a known bug reported in https://github.com/glue-viz/glue-astronomy/issues/52 ,
-aperture photometry and radial profile will report inaccurate results when you
-calculate them on dithered images linked by WCS *unless* you are on the reference image
-(this is usually the first loaded image).
-
 .. _known_issues_specviz:
 
 Specviz

jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.py

@@ -245,6 +245,8 @@ def vue_do_aper_phot(self, *args, **kwargs):
 
         data = self._selected_data
         reg = self._selected_subset
+        xcenter = reg.center.x
+        ycenter = reg.center.y
 
         # Reset last fitted model
         fit_model = None
@@ -257,6 +259,10 @@ def vue_do_aper_phot(self, *args, **kwargs):
                 bg = float(self.background_value)
             except ValueError:  # Clearer error message
                 raise ValueError('Missing or invalid background value')
+            if data.coords is not None:
+                sky_center = data.coords.pixel_to_world(xcenter, ycenter)
+            else:
+                sky_center = None
             aperture = regions2aperture(reg)
             include_pixarea_fac = False
             include_counts_fac = False
@@ -289,12 +295,10 @@ def vue_do_aper_phot(self, *args, **kwargs):
                     include_flux_scale = True
             phot_aperstats = ApertureStats(comp_data, aperture, wcs=data.coords, local_bkg=bg)
             phot_table = phot_aperstats.to_table(columns=(
-                'id', 'xcentroid', 'ycentroid', 'sky_centroid', 'sum', 'sum_aper_area',
+                'id', 'sum', 'sum_aper_area',
                 'min', 'max', 'mean', 'median', 'mode', 'std', 'mad_std', 'var',
                 'biweight_location', 'biweight_midvariance', 'fwhm', 'semimajor_sigma',
                 'semiminor_sigma', 'orientation', 'eccentricity'))  # Some cols excluded, add back as needed.  # noqa
-            phot_table['xcentroid'].unit = u.pix  # photutils only assumes, we make it real
-            phot_table['ycentroid'].unit = u.pix
             rawsum = phot_table['sum'][0]
 
             if include_pixarea_fac:
@@ -323,12 +327,14 @@ def vue_do_aper_phot(self, *args, **kwargs):
 
             # Extra info beyond photutils.
             phot_table.add_columns(
-                [bg, pixarea_fac, sum_ct, sum_ct_err, ctfac, sum_mag, flux_scale, data.label,
+                [xcenter * u.pix, ycenter * u.pix, sky_center,
+                 bg, pixarea_fac, sum_ct, sum_ct_err, ctfac, sum_mag, flux_scale, data.label,
                  reg.meta.get('label', ''), Time(datetime.utcnow())],
-                names=['background', 'pixarea_tot', 'aperture_sum_counts',
-                       'aperture_sum_counts_err', 'counts_fac', 'aperture_sum_mag', 'flux_scaling',
+                names=['xcenter', 'ycenter', 'sky_center', 'background', 'pixarea_tot',
+                       'aperture_sum_counts', 'aperture_sum_counts_err', 'counts_fac',
+                       'aperture_sum_mag', 'flux_scaling',
                        'data_label', 'subset_label', 'timestamp'],
-                indexes=[4, 6, 6, 6, 6, 6, 6, 21, 21, 21])
+                indexes=[1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 18, 18, 18])
 
             # Attach to app for Python extraction.
             if (not hasattr(self.app, '_aper_phot_results') or
@@ -352,9 +358,9 @@ def vue_do_aper_phot(self, *args, **kwargs):
             line_y_sc = bqplot.LinearScale()
 
             if self.current_plot_type == "Curve of Growth":
-                self._fig.title = 'Curve of growth from source centroid'
+                self._fig.title = 'Curve of growth from aperture center'
                 x_arr, sum_arr, x_label, y_label = _curve_of_growth(
-                    comp_data, phot_aperstats.centroid, aperture, phot_table['sum'][0],
+                    comp_data, (xcenter, ycenter), aperture, phot_table['sum'][0],
                     wcs=data.coords, background=bg, pixarea_fac=pixarea_fac)
                 self._fig.axes = [bqplot.Axis(scale=line_x_sc, label=x_label),
                                   bqplot.Axis(scale=line_y_sc, orientation='vertical',
@@ -370,33 +376,33 @@ def vue_do_aper_phot(self, *args, **kwargs):
                                               label=comp.units or 'Value')]
 
                 if self.current_plot_type == "Radial Profile":
-                    self._fig.title = 'Radial profile from source centroid'
+                    self._fig.title = 'Radial profile from aperture center'
                     x_data, y_data = _radial_profile(
-                        phot_aperstats.data_cutout, phot_aperstats.bbox, phot_aperstats.centroid,
+                        phot_aperstats.data_cutout, phot_aperstats.bbox, (xcenter, ycenter),
                         raw=False)
                     bqplot_line = bqplot.Lines(x=x_data, y=y_data, marker='circle',
                                                scales={'x': line_x_sc, 'y': line_y_sc},
                                                marker_size=32, colors='gray')
                 else:  # Radial Profile (Raw)
-                    self._fig.title = 'Raw radial profile from source centroid'
+                    self._fig.title = 'Raw radial profile from aperture center'
                     x_data, y_data = _radial_profile(
-                        phot_aperstats.data_cutout, phot_aperstats.bbox, phot_aperstats.centroid,
+                        phot_aperstats.data_cutout, phot_aperstats.bbox, (xcenter, ycenter),
                         raw=True)
                     bqplot_line = bqplot.Scatter(x=x_data, y=y_data, marker='circle',
                                                  scales={'x': line_x_sc, 'y': line_y_sc},
                                                  default_size=1, colors='gray')
 
                 # Fit Gaussian1D to radial profile data.
-                # mean is fixed at 0 because we recentered to centroid.
                 if self.fit_radial_profile:
                     fitter = LevMarLSQFitter()
                     y_max = y_data.max()
+                    x_mean = x_data[np.where(y_data == y_max)].mean()
                     std = 0.5 * (phot_table['semimajor_sigma'][0] +
                                  phot_table['semiminor_sigma'][0])
                     if isinstance(std, u.Quantity):
                         std = std.value
-                    gs = Gaussian1D(amplitude=y_max, mean=0, stddev=std,
-                                    fixed={'mean': True, 'amplitude': True},
+                    gs = Gaussian1D(amplitude=y_max, mean=x_mean, stddev=std,
+                                    fixed={'amplitude': True},
                                     bounds={'amplitude': (y_max * 0.5, y_max)})
                     if Version(astropy.__version__) < Version('5.2'):
                         fitter_kw = {}
@@ -433,13 +439,13 @@ def vue_do_aper_phot(self, *args, **kwargs):
                     continue
                 x = phot_table[key][0]
                 if (isinstance(x, (int, float, u.Quantity)) and
-                        key not in ('xcentroid', 'ycentroid', 'sky_centroid', 'sum_aper_area',
+                        key not in ('xcenter', 'ycenter', 'sky_center', 'sum_aper_area',
                                     'aperture_sum_counts', 'aperture_sum_mag')):
                     tmp.append({'function': key, 'result': f'{x:.4e}'})
-                elif key == 'sky_centroid' and x is not None:
-                    tmp.append({'function': 'RA centroid', 'result': f'{x.ra.deg:.6f} deg'})
-                    tmp.append({'function': 'Dec centroid', 'result': f'{x.dec.deg:.6f} deg'})
-                elif key in ('xcentroid', 'ycentroid', 'sum_aper_area'):
+                elif key == 'sky_center' and x is not None:
+                    tmp.append({'function': 'RA center', 'result': f'{x.ra.deg:.6f} deg'})
+                    tmp.append({'function': 'Dec center', 'result': f'{x.dec.deg:.6f} deg'})
+                elif key in ('xcenter', 'ycenter', 'sum_aper_area'):
                     tmp.append({'function': key, 'result': f'{x:.1f}'})
                 elif key == 'aperture_sum_counts' and x is not None:
                     tmp.append({'function': key, 'result':
@@ -452,7 +458,7 @@ def vue_do_aper_phot(self, *args, **kwargs):
             # Also display fit results
             fit_tmp = []
             if fit_model is not None and isinstance(fit_model, Gaussian1D):
-                for param in ('fwhm', 'amplitude'):  # mean is fixed at 0
+                for param in ('mean', 'fwhm', 'amplitude'):
                     p_val = getattr(fit_model, param)
                     if isinstance(p_val, Parameter):
                         p_val = p_val.value

jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.vue

@@ -141,7 +141,7 @@
       <jupyter-widget :widget="radial_plot" style="width: 100%; height: 480px" />
     </v-row>
 
-    <div v-if="plot_available && fit_radial_profile">
+    <div v-if="plot_available && fit_radial_profile && current_plot_type != 'Curve of Growth'">
       <j-plugin-section-header>Gaussian Fit Results</j-plugin-section-header>
       <v-row no-gutters>
         <v-col cols=6><U>Result</U></v-col>

jdaviz/configs/imviz/tests/test_parser.py

@@ -277,11 +277,11 @@ def test_parse_jwst_nircam_level2(self, imviz_helper):
         assert_allclose(phot_plugin.flux_scaling, 0.003631)
         phot_plugin.vue_do_aper_phot()
         tbl = imviz_helper.get_aperture_photometry_results()
-        assert_quantity_allclose(tbl[0]['xcentroid'], 970.935492 * u.pix)
-        assert_quantity_allclose(tbl[0]['ycentroid'], 1116.967619 * u.pix)
-        sky = tbl[0]['sky_centroid']
+        assert_quantity_allclose(tbl[0]['xcenter'], 970.95 * u.pix)
+        assert_quantity_allclose(tbl[0]['ycenter'], 1116.05 * u.pix)
+        sky = tbl[0]['sky_center']
         assert_allclose(sky.ra.deg, 80.48419863)
-        assert_allclose(sky.dec.deg, -69.494592)
+        assert_allclose(sky.dec.deg, -69.494608)
         data_unit = u.MJy / u.sr
         assert_quantity_allclose(tbl[0]['background'], 0.1741226315498352 * data_unit)
         assert_quantity_allclose(tbl[0]['sum'], 4.486487e-11 * u.MJy, rtol=1e-6)
@@ -402,9 +402,9 @@ def test_parse_hst_drz(self, imviz_helper):
         assert_allclose(phot_plugin.pixel_area, 0.0025)  # Not used but still auto-populated
         phot_plugin.vue_do_aper_phot()
         tbl = imviz_helper.get_aperture_photometry_results()
-        assert_quantity_allclose(tbl[0]['xcentroid'], 1487.60825422 * u.pix, atol=2 * u.pix)
-        assert_quantity_allclose(tbl[0]['ycentroid'], 2573.83983184 * u.pix, atol=2 * u.pix)
-        sky = tbl[0]['sky_centroid']
+        assert_quantity_allclose(tbl[0]['xcenter'], 1488.5 * u.pix, atol=2 * u.pix)
+        assert_quantity_allclose(tbl[0]['ycenter'], 2576 * u.pix, atol=2 * u.pix)
+        sky = tbl[0]['sky_center']
         assert_allclose(sky.ra.deg, 3.684062989070131, rtol=1e-3)
         assert_allclose(sky.dec.deg, 10.802045612042956, rtol=1e-3)
         data_unit = u.electron / u.s

jdaviz/configs/imviz/tests/test_simple_aper_phot.py

@@ -65,7 +65,7 @@ def test_plugin_wcs_dithered(self):
 
         # Check photometry results.
         assert tbl.colnames == [
-            'id', 'xcentroid', 'ycentroid', 'sky_centroid', 'background', 'sum',
+            'id', 'xcenter', 'ycenter', 'sky_center', 'background', 'sum',
             'sum_aper_area', 'pixarea_tot', 'aperture_sum_counts', 'aperture_sum_counts_err',
             'counts_fac', 'aperture_sum_mag', 'flux_scaling', 'min', 'max', 'mean', 'median',
             'mode', 'std', 'mad_std', 'var', 'biweight_location', 'biweight_midvariance',
@@ -99,11 +99,11 @@ def test_plugin_wcs_dithered(self):
         assert_array_equal(tbl['subset_label'], ['Subset 1', 'Subset 1'])
         assert tbl['timestamp'].scale == 'utc'
 
-        # Sky is the same but xcentroid different due to dithering.
+        # Sky is the same but xcenter different due to dithering.
         # The aperture sum is different too because mask is a little off limit in second image.
-        assert_quantity_allclose(tbl['xcentroid'], [4.5, 5.5] * u.pix)
-        assert_quantity_allclose(tbl['ycentroid'], 4.5 * u.pix)
-        sky = tbl['sky_centroid']
+        assert_quantity_allclose(tbl['xcenter'], [4.5, 5.5] * u.pix)
+        assert_quantity_allclose(tbl['ycenter'], 4.5 * u.pix)
+        sky = tbl['sky_center']
         assert_allclose(sky.ra.deg, 337.518943)
         assert_allclose(sky.dec.deg, -20.832083)
         assert_allclose(tbl['sum'], [63.617251, 62.22684693104279])
@@ -117,9 +117,9 @@ def test_plugin_wcs_dithered(self):
         tbl = self.imviz.get_aperture_photometry_results()
         assert len(tbl) == 3  # New result is appended
         assert tbl[-1]['id'] == 3
-        assert_quantity_allclose(tbl[-1]['xcentroid'], 4.5 * u.pix)
-        assert_quantity_allclose(tbl[-1]['ycentroid'], 2 * u.pix)
-        sky = tbl[-1]['sky_centroid']
+        assert_quantity_allclose(tbl[-1]['xcenter'], 4.5 * u.pix)
+        assert_quantity_allclose(tbl[-1]['ycenter'], 2 * u.pix)
+        sky = tbl[-1]['sky_center']
         assert_allclose(sky.ra.deg, 337.51894336144454)
         assert_allclose(sky.dec.deg, -20.832777499255897)
         assert_quantity_allclose(tbl[-1]['sum_aper_area'], 28.274334 * (u.pix * u.pix))
@@ -139,11 +139,11 @@ def test_plugin_wcs_dithered(self):
         tbl = self.imviz.get_aperture_photometry_results()
         assert len(tbl) == 4  # New result is appended
         assert tbl[-1]['id'] == 4
-        assert np.isnan(tbl[-1]['xcentroid'])
-        assert np.isnan(tbl[-1]['ycentroid'])
-        sky = tbl[-1]['sky_centroid']
-        assert np.isnan(sky.ra.deg)
-        assert np.isnan(sky.dec.deg)
+        assert_quantity_allclose(tbl[-1]['xcenter'], 4.5 * u.pix)
+        assert_quantity_allclose(tbl[-1]['ycenter'], 4.5 * u.pix)
+        sky = tbl[-1]['sky_center']
+        assert_allclose(sky.ra.deg, 337.51894336144454)
+        assert_allclose(sky.dec.deg, -20.832083)
         assert_quantity_allclose(tbl[-1]['sum_aper_area'], 81 * (u.pix * u.pix))
         assert_allclose(tbl[-1]['sum'], 0)
         assert_allclose(tbl[-1]['mean'], 0)
@@ -168,7 +168,7 @@ def test_plugin_wcs_dithered(self):
         # Curve of growth
         phot_plugin.current_plot_type = 'Curve of Growth'
         phot_plugin.vue_do_aper_phot()
-        assert phot_plugin._fig.title == 'Curve of growth from source centroid'
+        assert phot_plugin._fig.title == 'Curve of growth from aperture center'
 
 
 class TestSimpleAperPhot_NoWCS(BaseImviz_WCS_NoWCS):
@@ -187,7 +187,7 @@ def test_plugin_no_wcs(self):
         phot_plugin.vue_do_aper_phot()
         tbl = self.imviz.get_aperture_photometry_results()
         assert len(tbl) == 1  # Old table discarded due to incompatible column
-        assert_array_equal(tbl['sky_centroid'], None)
+        assert_array_equal(tbl['sky_center'], None)
 
 
 def test_annulus_background(imviz_helper):
@@ -263,6 +263,8 @@ def test_annulus_background(imviz_helper):
 
 # NOTE: Extracting the cutout for radial profile is aperture
 #       shape agnostic, so we use ellipse as representative case.
+# NOTE: This test only tests the radial profile algorithm and does
+#       not care if the actual plugin use centroid or not.
 class TestRadialProfile():
     def setup_class(self):
         data = np.ones((51, 51)) * u.nJy
@@ -286,6 +288,8 @@ def test_profile_imexam(self):
         assert_allclose(y_arr, 1)
 
 
+# NOTE: This test only tests the curve of growth algorithm and does
+#       not care if the actual plugin use centroid or not.
 @pytest.mark.parametrize('with_unit', (False, True))
 def test_curve_of_growth(with_unit):
     data = np.ones((51, 51))