Generalize the RHE Function beyond AIA

changelog/205.doc.rst

@@ -1,2 +1,2 @@
 Added a new example for `sunkit_image.radial.rhef` at :ref:`sphx_glr_generated_gallery_radial_histogram_equalization.py`
-Added RFEF to :ref:`sphx_glr_generated_gallery_radial_gradient_filters.py` as a comparison to the other filters.
+Added RHEF to :ref:`sphx_glr_generated_gallery_radial_gradient_filters.py` as a comparison to the other filters.

examples/radial_histogram_equalization.py

@@ -27,7 +27,7 @@
 radial_bin_edges = equally_spaced_bins(0, 2, aia_map.data.shape[0] // 2)
 radial_bin_edges *= u.R_sun
 
-rhef_map = radial.rhef(aia_map, radial_bin_edges, progress=False)
+rhef_map = radial.rhef(aia_map, radial_bin_edges=radial_bin_edges, progress=False)
 
 fig = plt.figure(figsize=(15, 10))
 

sunkit_image/radial.py

@@ -11,6 +11,7 @@
 import sunpy.map
 from sunpy.coordinates import frames
 
+
 from sunkit_image.utils import (
     apply_upsilon,
     bin_edge_summary,
@@ -104,7 +105,7 @@
 
 def intensity_enhance(
     smap,
-    radial_bin_edges,
+    radial_bin_edges=None,
     scale=None,
     summarize_bin_edges="center",
     summary=np.mean,
@@ -169,8 +170,8 @@
     `sunpy.map.Map`
         A SunPy map that has the emission above the normalization radius enhanced.
     """
-    # Get the radii for every pixel
-    map_r = find_pixel_radii(smap).to(u.R_sun)
+    # Handle the bin edges and radius array
+    radial_bin_edges, map_r = find_radial_bin_edges(smap, radial_bin_edges)
 
     # Get the radial intensity distribution
     radial_intensity = get_radial_intensity_summary(
@@ -207,21 +208,24 @@
     enhancement[map_r < normalization_radius] = 1
 
     # Return a map with the intensity enhanced above the normalization radius
+    # and the same meta data as the input map.
+
     new_map = sunpy.map.Map(smap.data * enhancement, smap.meta)
     new_map.plot_settings["norm"] = None
     return new_map
 
 
 def nrgf(
     smap,
-    radial_bin_edges,
+    radial_bin_edges=None,
     scale=None,
     intensity_summary=np.nanmean,
     intensity_summary_kwargs=None,
     width_function=np.std,
     width_function_kwargs=None,
     application_radius=1 * u.R_sun,
     progress=True,
+    fill=np.nan,
 ):
     """
     Implementation of the normalizing radial gradient filter (NRGF).
@@ -265,9 +269,11 @@
     application_radius : `astropy.units.Quantity`, optional
         The NRGF is applied to emission at radii above the application_radius.
         Defaults to 1 solar radii.
-    progress : `bool`, optional
-        Display a progressbar on the main loop.
-        Defaults to True.
+    progress : ``bool``, optional
+        Show a progressbar while computing
+    fill : ``any``, optional
+        The value to be placed outside of the bounds of the algorithm
+        Defaults to NAN.
 
     Returns
     -------
@@ -285,15 +291,9 @@
         width_function_kwargs = {}
     if intensity_summary_kwargs is None:
         intensity_summary_kwargs = {}
-    map_r = find_pixel_radii(smap).to(u.R_sun)
 
-    # To make sure bins are in the map.
-    if radial_bin_edges[1, -1] > np.max(map_r):
-        radial_bin_edges = equally_spaced_bins(
-            inner_value=radial_bin_edges[0, 0],
-            outer_value=np.max(map_r),
-            nbins=radial_bin_edges.shape[1],
-        )
+    # Handle the bin edges and radius array
+    radial_bin_edges, map_r = find_radial_bin_edges(smap, radial_bin_edges)
 
     # Radial intensity
     radial_intensity = get_radial_intensity_summary(
@@ -314,7 +314,7 @@
     )
 
     # Storage for the filtered data
-    data = np.zeros_like(smap.data)
+    data = np.ones_like(smap.data) * fill
 
     # Calculate the filter value for each radial bin.
     for i in tqdm(range(radial_bin_edges.shape[1]), desc="NRGF: ", disable=not progress):
@@ -403,6 +403,7 @@
     application_radius=1 * u.R_sun,
     number_angular_segments=130,
     progress=True,
+    fill=np.nan,
 ):
     """
     Implementation of the fourier normalizing radial gradient filter (FNRGF).
@@ -450,9 +451,11 @@
     number_angular_segments : `int`
         Number of angular segments in a circular annulus.
         Defaults to 130.
-    progress : `bool`, optional
-        Display a progressbar on the main loop.
-        Defaults to True.
+    progress : ``bool``, optional
+        Show a progressbar while computing
+    fill : ``any``, optional
+        The value to be placed outside of the bounds of the algorithm
+        Defaults to NAN.
 
     Returns
     -------
@@ -474,16 +477,8 @@
         msg = "Minimum value of order is 1"
         raise ValueError(msg)
 
-    # Get the radii for every pixel
-    map_r = find_pixel_radii(smap).to(u.R_sun)
-
-    # To make sure bins are in the map.
-    if radial_bin_edges[1, -1] > np.max(map_r):
-        radial_bin_edges = equally_spaced_bins(
-            inner_value=radial_bin_edges[0, 0],
-            outer_value=np.max(map_r),
-            nbins=radial_bin_edges.shape[1],
-        )
+    # Handle the bin edges and radius
+    radial_bin_edges, map_r = find_radial_bin_edges(smap, radial_bin_edges)
 
     # Get the Helioprojective coordinates of each pixel
     x, y = np.meshgrid(*[np.arange(v.value) for v in smap.dimensions]) * u.pix
@@ -499,7 +494,7 @@
     nbins = radial_bin_edges.shape[1]
 
     # Storage for the filtered data
-    data = np.zeros_like(smap.data)
+    data = np.ones_like(smap.data) * fill
 
     # Iterate over each circular ring
     for i in tqdm(range(nbins), desc="FNRGF: ", disable=not progress):
@@ -631,16 +626,37 @@
     return ranking_func
 
 
+def find_radial_bin_edges(smap, radial_bin_edges=None):
+    # Get the radii for every pixel, ensuring units are correct (in terms of pixels or solar radii)
+    map_r = find_pixel_radii(smap)
+    # Automatically generate radial bin edges if none are provided
+    if radial_bin_edges is None:
+        radial_bin_edges = equally_spaced_bins(0, np.max(map_r.value), smap.data.shape[0] // 2) * u.R_sun
+
+    # Ensure radial_bin_edges are within the bounds of the map_r values
+    if radial_bin_edges[1, -1] > np.max(map_r):
+        radial_bin_edges = (
+            equally_spaced_bins(
+                inner_value=radial_bin_edges[0, 0].to(u.R_sun).value,
+                outer_value=np.max(map_r.to(u.R_sun)).value,
+                nbins=radial_bin_edges.shape[1] // 2,
+            )
+            * u.R_sun
+        )
+    return radial_bin_edges, map_r
+
+
 @u.quantity_input(application_radius=u.R_sun, vignette=u.R_sun)
 def rhef(
     smap,
+    *,
     radial_bin_edges=None,
     application_radius=0 * u.R_sun,
     upsilon=0.35,
     method="numpy",
-    *,
-    vignette=1.5 * u.R_sun,
+    vignette=None,
     progress=True,
+    fill=np.nan,
 ):
     """
     Implementation of the Radial Histogram Equalizing Filter (RHEF).
@@ -658,33 +674,28 @@
     Parameters
     ----------
     smap : `sunpy.map.Map`
-        The sunpy map to enhance.
-    radial_bin_edges : `astropy.units.Quantity`
-        A two-dimensional array of bin edges of size ``[2, nbins]`` where ``nbins`` is
-        the number of bins.
+        The SunPy map to enhance using the RHEF algorithm.
+    radial_bin_edges : `astropy.units.Quantity`, optional
+        A two-dimensional array of bin edges of size ``[2, nbins]`` where ``nbins`` is the number of bins.
+        These define the radial segments where filtering is applied. If None, radial bins will be generated automatically.
     application_radius : `astropy.units.Quantity`, optional
-        The RHEF is applied to emission at radii above the application_radius.
+        The radius above which to apply the RHEF. Only regions with radii above this value will be filtered.
         Defaults to 0 solar radii.
-    upsilon : None, float, or tuple of `float`, optional
-        A double-sided gamma function applied to the equalized histograms.
-        See Equation (4.15) in the thesis.
-        Defaults to 0.35.
-    method : str
-    vignette: `astropy.units.Quantity`, optional
-        Set pixels above this radius to black.
-        Defaults to ``1.5*u.R_sun``.
-        If you want to disable this, pass in None.
-        Set pixels above this radius to black.
-        Defaults to None which is no vignette.
-        One suggested value is ``1.5*u.R_sun``.
-    progress: bool, optional
-        Display a progressbar on the main loop.
-        Defaults to True.
-
+    upsilon : float or None, optional
+        A double-sided gamma function to apply to modify the equalized histograms. Defaults to 0.35.
+    method : str, optional
+        Method used to rank the pixels for equalization. Defaults to 'inplace', with 'scipy' and 'numpy' as other options.
+    vignette : `astropy.units.Quantity`, optional
+        Radius beyond which pixels will be set to NAN. Defaults to None, must be in units that are compatible with "R_sun" as the value will be transformed.
+    progress : bool, optional
+        If True, display a progress bar during the filtering process. Defaults to True.
+    fill : ``any``, optional
+        The value to be placed outside of the bounds of the algorithm
+        Defaults to NAN.
     Returns
     -------
     `sunpy.map.Map`
-        A SunPy map that has had the RHEF applied to it.
+        A SunPy map with the Radial Histogram Equalizing Filter applied to it.
 
     References
     ----------
@@ -695,46 +706,36 @@
       https://www.proquest.com/docview/2759080511
     """
 
-    # Get the radii for every pixel
-    map_r = find_pixel_radii(smap).to(u.R_sun)
+    radial_bin_edges, map_r = find_radial_bin_edges(smap, radial_bin_edges)
 
-    if radial_bin_edges is None:
-        radial_bin_edges = equally_spaced_bins(0, 2, smap.data.shape[0] // 2)
-        radial_bin_edges *= u.R_sun
+    data = np.ones_like(smap.data) * fill
 
-    # Make sure bins are in the map.
-    if radial_bin_edges[1, -1] > np.max(map_r):
-        radial_bin_edges = equally_spaced_bins(
-            inner_value=radial_bin_edges[0, 0],
-            outer_value=np.max(map_r),
-            nbins=radial_bin_edges.shape[1],
-        )
-
-    # Allocate storage for the filtered data
-    data = np.zeros_like(smap.data)
-    meta = smap.meta
-
-    if radial_bin_edges.shape[1] > 2000:
-        progress = True
+    # Select the ranking method
+    ranking_func = _select_rank_method(method)
 
-    # Calculate the filter values for each radial bin.
+    # Loop over each radial bin to apply the filter
     for i in tqdm(range(radial_bin_edges.shape[1]), desc="RHEF: ", disable=not progress):
-        # Identify the appropriate radial slice
-        here = np.logical_and(map_r >= radial_bin_edges[0, i], map_r < radial_bin_edges[1, i])
+        # Identify pixels within the current radial bin
+        here = np.logical_and(
+            map_r >= radial_bin_edges[0, i].to(u.R_sun), map_r < radial_bin_edges[1, i].to(u.R_sun)
+        )
         if application_radius is not None and application_radius > 0:
             here = np.logical_and(here, map_r >= application_radius)
-
-        # Perform the filtering operation
-        ranking_func = _select_rank_method(method)
+        # Apply ranking function
         data[here] = ranking_func(smap.data[here])
         if upsilon is not None:
             data[here] = apply_upsilon(data[here], upsilon)
-    new_map = sunpy.map.Map(data, meta, autoalign=True)
+
+    new_map = sunpy.map.Map(data, smap.meta)
 
     if vignette is not None:
-        new_map = blackout_pixels_above_radius(new_map, vignette)
+        new_map = blackout_pixels_above_radius(new_map, vignette.to(u.R_sun))
+
+    # Adjust plot settings to remove extra normalization
+    # This must be done whenever one is adjusting
+    # the overall statistical distribution of values
 
-    # This must be done whenever one is adjusting the overall statistical distribution of values
     new_map.plot_settings["norm"] = None
 
+    # Return the new SunPy map with RHEF applied
     return new_map

sunkit_image/utils/utils.py

@@ -194,7 +194,9 @@
         msg = "Input array must be 2d!"
         raise ValueError(msg)
     if factor > 1:
-        congridx = RectBivariateSpline(np.arange(0, array.shape[0]), np.arange(0, array.shape[1]), array, kx=1, ky=1)
+        congridx = RectBivariateSpline(
+            np.arange(0, array.shape[0]), np.arange(0, array.shape[1]), array, kx=1, ky=1
+        )
         return congridx(np.arange(0, array.shape[0], 1 / factor), np.arange(0, array.shape[1], 1 / factor))
     return array
 
@@ -369,7 +371,7 @@
     return out_curve
 
 
-def blackout_pixels_above_radius(smap, radius_limit=1.5 * u.R_sun):
+def blackout_pixels_above_radius(smap, radius_limit=1.5 * u.R_sun, fill=np.nan):
     """
     Black out any pixels above a certain radius in a SunPy map.
 
@@ -379,6 +381,8 @@
         The input sunpy map.
     radius_limit : `astropy.units.Quantity`
         The radius limit above which to black out pixels.
+    fill : `any`
+        The value to use above the radius_limit.
 
     Returns
     -------
@@ -392,7 +396,7 @@
     mask = map_r > radius_limit
 
     # Apply the mask to the map data
-    masked_data = np.where(mask, 0, smap.data)
+    masked_data = np.where(mask, fill, smap.data)
 
     # Create a new map with the masked data
     return sunpy.map.Map(masked_data, smap.meta)