Merge pull request #1409 from pllim/fit-line-to-histogram

Simple Aperture Photometry: Gaussian1D fitting for radial profile

Author: pllim

CHANGES.rst

@@ -27,6 +27,10 @@ Cubeviz
 Imviz
 ^^^^^
 
+- Added the ability to fit Gaussian1D model to radial profile in
+  Simple Aperture Photometry plugin. Radial profile and curve of growth now center
+  on source centroid, not Subset center. [#1409]
+
 Mosviz
 ^^^^^^
 

docs/imviz/plugins.rst

@@ -160,7 +160,9 @@ an interactively selected region. A typical workflow is as follows:
     Caution: having too many data points may cause performance issues with this feature. 
     The exact limitations depend on your hardware.
 
-10. Once all inputs are populated correctly, click on the :guilabel:`CALCULATE`
+10. Toggle :guilabel:`Fit Gaussian` on to fit a `~astropy.modeling.functional_models.Gaussian1D`
+    model to the radial profile data. This is disabled for curve-of-growth.
+11. Once all inputs are populated correctly, click on the :guilabel:`CALCULATE`
     button to perform simple aperture photometry.
 
 .. note::
@@ -169,8 +171,8 @@ an interactively selected region. A typical workflow is as follows:
     However, if NaN exists in data, it will be treated as 0.
 
 When calculation is complete, a plot would show the radial profile
-of the background subtracted data and the photometry results are displayed under the
-:guilabel:`CALCULATE` button.
+of the background subtracted data and the photometry and model fitting (if requested)
+results are displayed under the :guilabel:`CALCULATE` button.
 
 .. figure:: img/imviz_radial_profile.png
     :alt: Imviz radial profile plot.
@@ -182,7 +184,14 @@ of the background subtracted data and the photometry results are displayed under
 
     Radial profile (raw).
 
-You can also retrieve the results as `~astropy.table.QTable` as follows,
+If you opted to fit a `~astropy.modeling.functional_models.Gaussian1D`
+to the radial profile, the last fitted model parameters will be displayed
+under the radial profile plot. The model itself can be obtained by as follows.
+See :ref:`astropy:astropy-modeling` on how to manipulate the model::
+
+    my_gaussian1d = imviz.app.fitted_models['phot_radial_profile']
+
+You can also retrieve the photometry results as `~astropy.table.QTable` as follows,
 assuming ``imviz`` is the instance of your Imviz application::
 
     results = imviz.get_aperture_photometry_results()

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

@@ -1,8 +1,13 @@
+import os
+import warnings
 from datetime import datetime
 
 import bqplot
 import numpy as np
 from astropy import units as u
+from astropy.modeling.fitting import LevMarLSQFitter
+from astropy.modeling import Parameter
+from astropy.modeling.models import Gaussian1D
 from astropy.table import QTable
 from astropy.time import Time
 from ipywidgets import widget_serialization
@@ -41,6 +46,8 @@ class SimpleAperturePhotometry(TemplateMixin, DatasetSelectMixin):
     current_plot_type = Unicode().tag(sync=True)
     plot_available = Bool(False).tag(sync=True)
     radial_plot = Any('').tag(sync=True, **widget_serialization)
+    fit_radial_profile = Bool(False).tag(sync=True)
+    fit_results = List().tag(sync=True)
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
@@ -63,6 +70,7 @@ def __init__(self, *args, **kwargs):
         self._fig = bqplot.Figure()
         self.plot_types = ["Curve of Growth", "Radial Profile", "Radial Profile (Raw)"]
         self.current_plot_type = self.plot_types[0]
+        self._fitted_model_name = 'phot_radial_profile'
 
     def reset_results(self):
         self.result_available = False
@@ -218,6 +226,11 @@ def vue_do_aper_phot(self, *args, **kwargs):
         data = self._selected_data
         reg = self._selected_subset
 
+        # Reset last fitted model
+        fit_model = None
+        if self._fitted_model_name in self.app.fitted_models:
+            del self.app.fitted_models[self._fitted_model_name]
+
         try:
             comp = data.get_component(data.main_components[0])
             try:
@@ -319,38 +332,68 @@ 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 Subset center'
+                self._fig.title = 'Curve of growth from source centroid'
                 x_arr, sum_arr, x_label, y_label = _curve_of_growth(
-                    comp_data, aperture, phot_table['sum'][0], wcs=data.coords,
-                    background=bg, pixarea_fac=pixarea_fac)
+                    comp_data, phot_aperstats.centroid, 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',
                                               label=y_label)]
                 bqplot_line = bqplot.Lines(x=x_arr, y=sum_arr, marker='circle',
                                            scales={'x': line_x_sc, 'y': line_y_sc},
                                            marker_size=32, colors='gray')
+                bqplot_marks = [bqplot_line]
 
             else:  # Radial profile
                 self._fig.axes = [bqplot.Axis(scale=line_x_sc, label='pix'),
                                   bqplot.Axis(scale=line_y_sc, orientation='vertical',
                                               label=comp.units or 'Value')]
 
                 if self.current_plot_type == "Radial Profile":
-                    self._fig.title = 'Radial profile from Subset center'
+                    self._fig.title = 'Radial profile from source centroid'
                     x_data, y_data = _radial_profile(
-                        phot_aperstats.data_cutout, phot_aperstats.bbox, aperture, raw=False)
+                        phot_aperstats.data_cutout, phot_aperstats.bbox, phot_aperstats.centroid,
+                        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 Subset center'
-                    radial_r, radial_img = _radial_profile(
-                        phot_aperstats.data_cutout, phot_aperstats.bbox, aperture, raw=True)
-                    bqplot_line = bqplot.Scatter(x=radial_r, y=radial_img, marker='circle',
+                    self._fig.title = 'Raw radial profile from source centroid'
+                    x_data, y_data = _radial_profile(
+                        phot_aperstats.data_cutout, phot_aperstats.bbox, phot_aperstats.centroid,
+                        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')
 
-            self._fig.marks = [bqplot_line]
+                # 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()
+                    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},
+                                    bounds={'amplitude': (y_max * 0.5, y_max)})
+                    with warnings.catch_warnings(record=True) as warns:
+                        fit_model = fitter(gs, x_data, y_data)
+                    if len(warns) > 0:
+                        msg = os.linesep.join([str(w.message) for w in warns])
+                        self.hub.broadcast(SnackbarMessage(
+                            f"Radial profile fitting: {msg}", color='warning', sender=self))
+                    y_fit = fit_model(x_data)
+                    self.app.fitted_models[self._fitted_model_name] = fit_model
+                    bqplot_fit = bqplot.Lines(x=x_data, y=y_fit, marker=None,
+                                              scales={'x': line_x_sc, 'y': line_y_sc},
+                                              colors='magenta', line_style='dashed')
+                    bqplot_marks = [bqplot_line, bqplot_fit]
+                else:
+                    bqplot_marks = [bqplot_line]
+
+            self._fig.marks = bqplot_marks
 
         except Exception as e:  # pragma: no cover
             self.reset_results()
@@ -379,7 +422,18 @@ def vue_do_aper_phot(self, *args, **kwargs):
                                 f'{x:.4e} ({phot_table["aperture_sum_counts_err"][0]:.4e})'})
                 else:
                     tmp.append({'function': key, 'result': str(x)})
+
+            # 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
+                    p_val = getattr(fit_model, param)
+                    if isinstance(p_val, Parameter):
+                        p_val = p_val.value
+                    fit_tmp.append({'function': param, 'result': f'{p_val:.4e}'})
+
             self.results = tmp
+            self.fit_results = fit_tmp
             self.result_available = True
             self.radial_plot = self._fig
             self.bqplot_figs_resize = [self._fig]
@@ -389,7 +443,7 @@ def vue_do_aper_phot(self, *args, **kwargs):
 # NOTE: These are hidden because the APIs are for internal use only
 # but we need them as a separate functions for unit testing.
 
-def _radial_profile(radial_cutout, reg_bb, aperture, raw=False):
+def _radial_profile(radial_cutout, reg_bb, centroid, raw=False):
     """Calculate radial profile.
 
     Parameters
@@ -400,23 +454,24 @@ def _radial_profile(radial_cutout, reg_bb, aperture, raw=False):
     reg_bb : obj
         Bounding box from ``ApertureStats``.
 
-    aperture : obj
-        ``photutils`` aperture object.
+    centroid : tuple of int
+        ``ApertureStats`` centroid or desired center in ``(x, y)``.
 
     raw : bool
         If `True`, returns raw data points for scatter plot.
         Otherwise, use ``imexam`` algorithm for a clean plot.
 
     """
     reg_ogrid = np.ogrid[reg_bb.iymin:reg_bb.iymax, reg_bb.ixmin:reg_bb.ixmax]
-    radial_dx = reg_ogrid[1] - aperture.positions[0]
-    radial_dy = reg_ogrid[0] - aperture.positions[1]
+    radial_dx = reg_ogrid[1] - centroid[0]
+    radial_dy = reg_ogrid[0] - centroid[1]
     radial_r = np.hypot(radial_dx, radial_dy)[~radial_cutout.mask].ravel()  # pix
     radial_img = radial_cutout.compressed()  # data unit
 
     if raw:
-        x_arr = radial_r
-        y_arr = radial_img
+        i_arr = np.argsort(radial_r)
+        x_arr = radial_r[i_arr]
+        y_arr = radial_img[i_arr]
     else:
         # This algorithm is from the imexam package,
         # see licenses/IMEXAM_LICENSE.txt for more details
@@ -427,7 +482,7 @@ def _radial_profile(radial_cutout, reg_bb, aperture, raw=False):
     return x_arr, y_arr
 
 
-def _curve_of_growth(data, aperture, final_sum, wcs=None, background=0, n_datapoints=10,
+def _curve_of_growth(data, centroid, aperture, final_sum, wcs=None, background=0, n_datapoints=10,
                      pixarea_fac=None):
     """Calculate curve of growth for aperture photometry.
 
@@ -436,8 +491,14 @@ def _curve_of_growth(data, aperture, final_sum, wcs=None, background=0, n_datapo
     data : ndarray or `~astropy.units.Quantity`
         Data for the calculation.
 
+    centroid : tuple of int
+        ``ApertureStats`` centroid or desired center in ``(x, y)``.
+
     aperture : obj
-        ``photutils`` aperture object.
+        ``photutils`` aperture to use, except its center will be
+        changed to the given ``centroid``. This is because the aperture
+        might be hand-drawn and a more accurate centroid has been
+        recalculated separately.
 
     final_sum : float or `~astropy.units.Quantity`
         Aperture sum that is already calculated in the
@@ -477,20 +538,20 @@ def _curve_of_growth(data, aperture, final_sum, wcs=None, background=0, n_datapo
     if isinstance(aperture, CircularAperture):
         x_label = 'Radius (pix)'
         x_arr = np.linspace(0, aperture.r, num=n_datapoints)[1:]
-        aper_list = [CircularAperture(aperture.positions, cur_r) for cur_r in x_arr[:-1]]
+        aper_list = [CircularAperture(centroid, cur_r) for cur_r in x_arr[:-1]]
     elif isinstance(aperture, EllipticalAperture):
         x_label = 'Semimajor axis (pix)'
         x_arr = np.linspace(0, aperture.a, num=n_datapoints)[1:]
         a_arr = x_arr[:-1]
         b_arr = aperture.b * a_arr / aperture.a
-        aper_list = [EllipticalAperture(aperture.positions, cur_a, cur_b, theta=aperture.theta)
+        aper_list = [EllipticalAperture(centroid, cur_a, cur_b, theta=aperture.theta)
                      for (cur_a, cur_b) in zip(a_arr, b_arr)]
     elif isinstance(aperture, RectangularAperture):
         x_label = 'Width (pix)'
         x_arr = np.linspace(0, aperture.w, num=n_datapoints)[1:]
         w_arr = x_arr[:-1]
         h_arr = aperture.h * w_arr / aperture.w
-        aper_list = [RectangularAperture(aperture.positions, cur_w, cur_h, theta=aperture.theta)
+        aper_list = [RectangularAperture(centroid, cur_w, cur_h, theta=aperture.theta)
                      for (cur_w, cur_h) in zip(w_arr, h_arr)]
     else:
         raise TypeError(f'Unsupported aperture: {aperture}')

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

@@ -110,6 +110,16 @@
           ></v-select>
         </v-row>
 
+        <v-row v-if="current_plot_type.indexOf('Radial Profile') != -1">
+
+          <v-switch
+            label="Fit Gaussian"
+            hint="Fit Gaussian1D to radial profile"
+            v-model="fit_radial_profile"
+            persistent-hint>
+          </v-switch>
+        </v-row>
+
         <v-row justify="end">
           <v-btn color="primary" text @click="do_aper_phot">Calculate</v-btn>
         </v-row>
@@ -123,8 +133,25 @@
       <jupyter-widget :widget="radial_plot" style="width: 100%; height: 480px" />
     </v-row>
 
+    <div v-if="plot_available && fit_radial_profile">
+      <j-plugin-section-header>Gaussian Fit Results</j-plugin-section-header>
+      <v-row no-gutters>
+        <v-col cols=6><U>Result</U></v-col>
+        <v-col cols=6><U>Value</U></v-col>
+      </v-row>
+      <v-row
+        v-for="item in fit_results"
+        :key="item.function"
+        no-gutters>
+        <v-col cols=6>
+          {{  item.function  }}
+        </v-col>
+        <v-col cols=6>{{ item.result }}</v-col>
+      </v-row>
+    </div>
+
     <div v-if="result_available">
-      <j-plugin-section-header>Results</j-plugin-section-header>
+      <j-plugin-section-header>Photometry Results</j-plugin-section-header>
       <v-row no-gutters>
         <v-col cols=6><U>Result</U></v-col>
         <v-col cols=6><U>Value</U></v-col>

jdaviz/configs/imviz/tests/test_simple_aper_phot.py

@@ -18,6 +18,10 @@ def test_plugin_wcs_dithered(self):
 
         phot_plugin = self.imviz.app.get_tray_item_from_name('imviz-aper-phot-simple')
 
+        # Model fitting is already tested in astropy.
+        # Here, we enable it just to make sure it does not crash.
+        phot_plugin.fit_radial_profile = True
+
         # Make sure invalid Data/Subset selection does not crash plugin.
         with pytest.raises(ValueError):
             phot_plugin.dataset_selected = 'no_such_data'
@@ -164,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 Subset center'
+        assert phot_plugin._fig.title == 'Curve of growth from source centroid'
 
 
 class TestSimpleAperPhot_NoWCS(BaseImviz_WCS_NoWCS):
@@ -255,21 +259,22 @@ def test_annulus_background(imviz_helper):
 class TestRadialProfile():
     def setup_class(self):
         data = np.ones((51, 51)) * u.nJy
-        self.aperture = EllipticalAperture((25, 25), 20, 15)
-        phot_aperstats = ApertureStats(data, self.aperture)
+        aperture = EllipticalAperture((25, 25), 20, 15)
+        phot_aperstats = ApertureStats(data, aperture)
         self.data_cutout = phot_aperstats.data_cutout
         self.bbox = phot_aperstats.bbox
+        self.centroid = phot_aperstats.centroid
 
     def test_profile_raw(self):
-        x_arr, y_arr = _radial_profile(self.data_cutout, self.bbox, self.aperture, raw=True)
+        x_arr, y_arr = _radial_profile(self.data_cutout, self.bbox, self.centroid, raw=True)
         # Too many data points to compare each one for X.
         assert x_arr.shape == y_arr.shape == (923, )
         assert_allclose(x_arr.min(), 0)
         assert_allclose(x_arr.max(), 19.4164878389476)
         assert_allclose(y_arr, 1)
 
     def test_profile_imexam(self):
-        x_arr, y_arr = _radial_profile(self.data_cutout, self.bbox, self.aperture, raw=False)
+        x_arr, y_arr = _radial_profile(self.data_cutout, self.bbox, self.centroid, raw=False)
         assert_allclose(x_arr, range(20))
         assert_allclose(y_arr, 1)
 
@@ -293,11 +298,12 @@ def test_curve_of_growth(with_unit):
                  RectangularAperture(cen, 20, 15))
 
     for aperture in apertures:
-        final_sum = ApertureStats(data, aperture).sum
+        astat = ApertureStats(data, aperture)
+        final_sum = astat.sum
         if pixarea_fac is not None:
             final_sum = final_sum * pixarea_fac
         x_arr, sum_arr, x_label, y_label = _curve_of_growth(
-            data, aperture, final_sum, background=bg, pixarea_fac=pixarea_fac)
+            data, astat.centroid, aperture, final_sum, background=bg, pixarea_fac=pixarea_fac)
         assert_allclose(x_arr, [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
         assert y_label == expected_ylabel
 
@@ -316,5 +322,5 @@ def test_curve_of_growth(with_unit):
             assert_allclose(sum_arr, [3, 12, 27, 48, 75, 108, 147, 192, 243, 300])
 
     with pytest.raises(TypeError, match='Unsupported aperture'):
-        _curve_of_growth(data, EllipticalAnnulus(cen, 3, 8, 5), 100,
+        _curve_of_growth(data, cen, EllipticalAnnulus(cen, 3, 8, 5), 100,
                          pixarea_fac=pixarea_fac)