Fix center profile

stellarphot/photometry/profiles.py

@@ -178,8 +178,10 @@ def __init__(
             profile_radius = cutout_size // 2
 
         radii = np.linspace(0, profile_radius, profile_radius + 1)
-        # Get a rough profile without background subtraction
-        self._radial_profile = RadialProfile(self._data, self._cen, radii)
+        # Get a rough profile with rough background subtraction -- note that
+        # NO background subtraction does not work.
+        background = sigma_clipped_stats(self._cutout.data)[1]
+        self._radial_profile = RadialProfile(self._data - background, self._cen, radii)
 
         self._sky_area = None
 
@@ -245,17 +247,32 @@ def pixel_values_in_profile(self):
         """
         radii = []
         pixel_values = []
-        for rad, ap in zip(
-            self.radial_profile.radius, self.radial_profile.apertures, strict=True
-        ):
+        for ap in self.radial_profile.apertures:
+            # Calculate the distance of each pixel from the center
+            grid_x, grid_y = np.mgrid[: self._data.shape[0], : self._data.shape[1]]
+            dist_from_cen = np.sqrt(
+                (grid_x - self._cen[1]) ** 2 + (grid_y - self._cen[0]) ** 2
+            )
+
+            # Get only the data in this aperture
             ap_mask = ap.to_mask(method="center")
             ap_data = ap_mask.multiply(self._data)
+            dist_from_cen = ap_mask.multiply(dist_from_cen)
+
+            # Drop any data where the aperture mask was zero and flatten
             good_data = ap_data != 0
-            pixel_values.extend(ap_data[good_data].flatten())
-            radii.extend([rad] * good_data.sum())
+            ap_exact_radius = dist_from_cen[good_data].flatten()
+            ap_data = ap_data[good_data].flatten()
+
+            # Sort so that plots look ok
+            sorted_radii = np.argsort(ap_exact_radius)
+            radii.extend(ap_exact_radius[sorted_radii])
+            pixel_values.extend(ap_data[sorted_radii])
+
         radii = np.array(radii)
         pixel_values = np.array(pixel_values)
-        return radii, pixel_values
+        # Subtract the background
+        return radii, pixel_values - self.sky_pixel_value
 
     @lazyproperty
     def curve_of_growth(self):

stellarphot/photometry/tests/test_profiles.py

@@ -72,6 +72,8 @@ def test_find_center_no_star():
 
 
 def test_radial_profile():
+    # Test that both curve of growth and radial profile are correct
+
     image = make_gaussian_sources_image(SHAPE, STARS)
     for row in STARS:
         cen = find_center(image, (row["x_mean"], row["y_mean"]), max_iters=10)
@@ -81,13 +83,21 @@ def test_radial_profile():
         # cutout size around 60.
         rad_prof = CenterAndProfile(image, cen, cutout_size=60, profile_radius=30)
 
+        # Test that the curve of growth is correct
+
         # Numerical value below is integral of input 2D gaussian, 2pi A sigma^2
         expected_integral = 2 * np.pi * row["amplitude"] * row["x_stddev"] ** 2
-
         np.testing.assert_allclose(
             rad_prof.curve_of_growth.profile[-1], expected_integral, atol=50
         )
 
+        # Test that the radial profile is correct by comparing pixel values to a
+        # gaussian fit to the profile.
+        data_radii, data_counts = rad_prof.pixel_values_in_profile
+        expected_profile = rad_prof.radial_profile.gaussian_fit(data_radii)
+
+        np.testing.assert_allclose(data_counts, expected_profile, atol=20)
+
 
 def test_radial_profile_exposure_is_nan():
     # Check that using an exposure value of NaN returns NaN for the SNR and noise
@@ -102,3 +112,32 @@ def test_radial_profile_exposure_is_nan():
     assert np.isnan(rad_prof.noise(c, np.nan)[-1])
     assert all(np.isfinite(rad_prof.curve_of_growth.profile))
     assert all(np.isfinite(rad_prof.radial_profile.profile))
+
+
+def test_radial_profile_with_background():
+    # Regression test for #328 -- image with a background level
+    image = make_gaussian_sources_image(SHAPE, STARS)
+    image = image + make_noise_image(
+        image.shape, distribution="gaussian", mean=100, stddev=0.1
+    )
+    for row in STARS:
+        cen = find_center(image, (row["x_mean"], row["y_mean"]), max_iters=10)
+
+        # The "stars" have FWHM around 9.5, so make the cutouts used for finding the
+        # stars fairly big -- the bare minimum would be a radius of 3 FWHM, which is a
+        # cutout size around 60.
+        rad_prof = CenterAndProfile(image, cen, cutout_size=60, profile_radius=30)
+
+        # Numerical value below is integral of input 2D gaussian, 2pi A sigma^2
+        expected_integral = 2 * np.pi * row["amplitude"] * row["x_stddev"] ** 2
+
+        np.testing.assert_allclose(
+            rad_prof.curve_of_growth.profile[-1], expected_integral, atol=50
+        )
+
+        # Test that the radial profile is correct by comparing pixel values to a
+        # gaussian fit to the profile.
+        data_radii, data_counts = rad_prof.pixel_values_in_profile
+        expected_profile = rad_prof.radial_profile.gaussian_fit(data_radii)
+
+        np.testing.assert_allclose(data_counts, expected_profile, atol=20)

stellarphot/settings/custom_widgets.py

@@ -14,6 +14,8 @@
     ui_generator,
 )
 
+__all__ = ["ChooseOrMakeNew", "Confirm"]
+
 DEFAULT_BUTTON_WIDTH = "300px"
 
 
@@ -39,7 +41,9 @@ class ChooseOrMakeNew(ipw.VBox):
         PassbandMap.__name__,
     ]
 
-    def __init__(self, item_type_name, *arg, _testing_path=None, **kwargs):
+    def __init__(
+        self, item_type_name, *arg, details_hideable=False, _testing_path=None, **kwargs
+    ):
         if item_type_name not in self._known_types:
             raise ValueError(
                 f"Unknown item type {item_type_name}. Must "
@@ -60,15 +64,26 @@ def __init__(self, item_type_name, *arg, _testing_path=None, **kwargs):
         # and track if we are making a new item
         self._making_new = False
 
+        # keep track of whether there is a "show details" checkbox
+        self._show_details_shown = details_hideable
+
         self._display_name = item_type_name.replace("_", " ")
 
         # Create the child widgets
+
+        # Descriptive title
         self._title = ipw.HTML(
             value=(f"<h2>Choose a {self._display_name} " "or make a new one</h2>")
         )
 
+        # Selector for existing items or to make a new one
         self._choose_existing = ipw.Dropdown(description="")
 
+        # Option to show/hide details, only displayed if user wants it.
+        self._show_details_ui = ipw.Checkbox(description="Show details", value=True)
+        self._show_details_ui.layout.display = "flex" if details_hideable else "none"
+        self._show_details_cached_value = self._show_details_ui.value
+
         self._edit_delete_container = ipw.HBox(
             # width below was chosen to match the dropdown...would prefer to
             # determine this programmatically but don't know how.
@@ -83,19 +98,28 @@ def __init__(self, item_type_name, *arg, _testing_path=None, **kwargs):
             description=f"Delete this {self._display_name}",
         )
 
+        # Put almost everything into a VBox
+        self._details_box = ipw.VBox()
+
         self._edit_delete_container.children = [self._edit_button, self._delete_button]
 
         self._confirm_edit_delete = Confirm()
 
         self._item_widget, self._widget_value_new_item = self._make_new_widget()
 
+        # Put all of the details into a box that can be easily hidden
+        self._details_box.children = [
+            self._edit_delete_container,
+            self._confirm_edit_delete,
+            self._item_widget,
+        ]
+
         # Build the main widget
         self.children = [
             self._title,
             self._choose_existing,
-            self._edit_delete_container,
-            self._confirm_edit_delete,
-            self._item_widget,
+            self._show_details_ui,
+            self._details_box,
         ]
 
         # Set up the dropdown widget
@@ -138,6 +162,17 @@ def __init__(self, item_type_name, *arg, _testing_path=None, **kwargs):
         # Respond when user wants to make a new thing
         self._choose_existing.observe(self._handle_selection, names="value")
 
+        # Set up an observer to show/hide the details box if the check box
+        # is clicked
+        self._show_details_ui.observe(self._show_details_handler, names="value")
+
+    @property
+    def value(self):
+        """
+        The value of the widget.
+        """
+        return self._item_widget.model(**self._item_widget.value)
+
     def _save_confirmation(self):
         """
         Function to attach to the save button to show the confirmation widget if
@@ -173,6 +208,12 @@ def _handle_selection(self, change):
             # Hide the edit button
             self._edit_delete_container.layout.display = "none"
 
+            # Make sure details are shown and hide the "show details" checkbox
+            if self._show_details_shown:
+                self._show_details_cached_value = self._show_details_ui.value
+                self._show_details_ui.value = True
+                self._show_details_ui.layout.display = "none"
+
             # This sets the ui back to its original state when created, i.e.
             # everything is empty.
             self._item_widget._init_ui()
@@ -217,6 +258,13 @@ def _handle_selection(self, change):
             # but does no harm in the other cases
             self._set_disable_state_nested_models(self._item_widget, True)
 
+            # We may have arrived here by choosing a different item while
+            # making a new one, so we restore the state of the "show details"
+            # checkbox.
+            if self._show_details_shown:
+                self._show_details_ui.layout.display = "flex"
+                self._show_details_ui.value = self._show_details_cached_value
+
     def _edit_button_action(self, _):
         """
         Handle the edit button being clicked.
@@ -264,6 +312,16 @@ def _delete_button_action(self, _):
         self._set_confirm_message()
         self._confirm_edit_delete.show()
 
+    def _show_details_handler(self, change):
+        """
+        Show or hide the details box based on the value of the checkbox.
+        """
+        if self._show_details_ui.layout.display == "none":
+            # The element is hidden, so just return
+            return
+
+        self._details_box.layout.display = "flex" if change["new"] else "none"
+
     def _set_disable_state_nested_models(self, top, value):
         """
         When a one model contains another and the top-level model widget
@@ -338,6 +396,9 @@ def saver():
             else:
                 # If saving works, we update the choices and select the new item
                 self._making_new = False
+                if self._show_details_shown:
+                    self._show_details_ui.layout.display = "flex"
+                    self._show_details_ui.value = self._show_details_cached_value
                 update_choices_and_select_new()
 
         def update_choices_and_select_new():
@@ -424,6 +485,13 @@ def confirm_handler(change):
                             self._choose_existing.value = self._choose_existing.options[
                                 0
                             ][1]
+                    if self._making_new:
+                        if self._show_details_shown:
+                            self._show_details_ui.layout.display = "flex"
+                            self._show_details_ui.value = (
+                                self._show_details_cached_value
+                            )
+
                     # We are done editing/making new regardless of
                     # the confirmation outcome
                     self._editing = False