From cb010a993828adfc5e7e19af82cccf12acf10a6a Mon Sep 17 00:00:00 2001 From: Nabil Freij Date: Tue, 29 Oct 2024 13:48:45 -0700 Subject: [PATCH] Update code --- changelog/205.doc.rst | 1 + changelog/231.bugfix.rst | 1 + sunkit_image/radial.py | 134 ++++++++++++++++-------------------- sunkit_image/utils/utils.py | 43 ++++++++++++ 4 files changed, 103 insertions(+), 76 deletions(-) create mode 100644 changelog/231.bugfix.rst diff --git a/changelog/205.doc.rst b/changelog/205.doc.rst index 3e85ba65..0a4ed021 100644 --- a/changelog/205.doc.rst +++ b/changelog/205.doc.rst @@ -1,2 +1,3 @@ Added a new example for `sunkit_image.radial.rhef` at :ref:`sphx_glr_generated_gallery_radial_histogram_equalization.py` + Added RHEF to :ref:`sphx_glr_generated_gallery_radial_gradient_filters.py` as a comparison to the other filters. diff --git a/changelog/231.bugfix.rst b/changelog/231.bugfix.rst new file mode 100644 index 00000000..887e5583 --- /dev/null +++ b/changelog/231.bugfix.rst @@ -0,0 +1 @@ +Improved the Radial Histogram Equalizing Filter (RHEF) (`sunkit_image.radial.rhef`) to work with images outside of SDO/AIA. diff --git a/sunkit_image/radial.py b/sunkit_image/radial.py index 1cd3103d..fca007fb 100644 --- a/sunkit_image/radial.py +++ b/sunkit_image/radial.py @@ -11,13 +11,11 @@ import sunpy.map from sunpy.coordinates import frames - from sunkit_image.utils import ( apply_upsilon, bin_edge_summary, blackout_pixels_above_radius, - equally_spaced_bins, - find_pixel_radii, + find_radial_bin_edges, get_radial_intensity_summary, ) @@ -102,9 +100,39 @@ def _normalize_fit_radial_intensity(radii, polynomial, normalization_radius): polynomial, ) +def _select_rank_method(method): + # For now, we have more than one option for ranking the values + def _percentile_ranks_scipy(arr): + from scipy import stats + + return stats.rankdata(arr, method="average") / len(arr) + + def _percentile_ranks_numpy(arr): + sorted_indices = np.argsort(arr) + ranks = np.empty_like(sorted_indices) + ranks[sorted_indices] = np.arange(1, len(arr) + 1) + return ranks / float(len(arr)) + + def _percentile_ranks_numpy_inplace(arr): + sorted_indices = np.argsort(arr) + arr[sorted_indices] = np.arange(1, len(arr) + 1) + return arr / float(len(arr)) + + # Select the sort method + if method == "inplace": + ranking_func = _percentile_ranks_numpy_inplace + elif method == "numpy": + ranking_func = _percentile_ranks_numpy + elif method == "scipy": + ranking_func = _percentile_ranks_scipy + else: + msg = f"{method} is invalid. Allowed values are 'inplace', 'numpy', 'scipy'" + raise NotImplementedError(msg) + return ranking_func def intensity_enhance( smap, + *, radial_bin_edges=None, scale=None, summarize_bin_edges="center", @@ -136,9 +164,10 @@ def intensity_enhance( ---------- smap : `sunpy.map.Map` The sunpy map to enhance. - radial_bin_edges : `astropy.units.Quantity` + 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. + Defaults to `None` which will use equally spaced bins. scale : `astropy.units.Quantity`, optional The radius of the Sun expressed in map units. For example, in typical Helioprojective Cartesian maps the solar radius is expressed in @@ -209,7 +238,6 @@ def intensity_enhance( # 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 @@ -224,7 +252,7 @@ def nrgf( width_function=np.std, width_function_kwargs=None, application_radius=1 * u.R_sun, - progress=True, + progress=False, fill=np.nan, ): """ @@ -269,11 +297,12 @@ def nrgf( 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 - Show a progressbar while computing - fill : ``any``, optional - The value to be placed outside of the bounds of the algorithm - Defaults to NAN. + progress : `bool`, optional + Show a progressbar while computing. + Defaults to `False`. + fill : Any, optional + The value to be placed outside of the bounds of the algorithm. + Defaults to NaN. Returns ------- @@ -402,7 +431,7 @@ def fnrgf( width_function=np.std, application_radius=1 * u.R_sun, number_angular_segments=130, - progress=True, + progress=False, fill=np.nan, ): """ @@ -451,11 +480,12 @@ def fnrgf( number_angular_segments : `int` Number of angular segments in a circular annulus. Defaults to 130. - 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. + progress : `bool`, optional + Show a progressbar while computing. + Defaults to `False`. + fill : Any, optional + The value to be placed outside of the bounds of the algorithm. + Defaults to NaN. Returns ------- @@ -595,57 +625,6 @@ def fnrgf( return new_map -def _select_rank_method(method): - # For now, we have more than one option for ranking the values - def _percentile_ranks_scipy(arr): - from scipy import stats - - return stats.rankdata(arr, method="average") / len(arr) - - def _percentile_ranks_numpy(arr): - sorted_indices = np.argsort(arr) - ranks = np.empty_like(sorted_indices) - ranks[sorted_indices] = np.arange(1, len(arr) + 1) - return ranks / float(len(arr)) - - def _percentile_ranks_numpy_inplace(arr): - sorted_indices = np.argsort(arr) - arr[sorted_indices] = np.arange(1, len(arr) + 1) - return arr / float(len(arr)) - - # Select the sort method - if method == "inplace": - ranking_func = _percentile_ranks_numpy_inplace - elif method == "numpy": - ranking_func = _percentile_ranks_numpy - elif method == "scipy": - ranking_func = _percentile_ranks_scipy - else: - msg = f"{method} is invalid. Allowed values are 'inplace', 'numpy', 'scipy'" - raise NotImplementedError(msg) - 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, @@ -655,7 +634,7 @@ def rhef( upsilon=0.35, method="numpy", vignette=None, - progress=True, + progress=False, fill=np.nan, ): """ @@ -683,15 +662,18 @@ def rhef( Defaults to 0 solar radii. 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. + method : {"inplace", "numpy", "scipy"}, optional + Method used to rank the pixels for equalization. + Defaults to 'inplace'. 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. + progress : `bool`, optional + Show a progressbar while computing. + Defaults to `False`. + fill : Any, optional + The value to be placed outside of the bounds of the algorithm. + Defaults to NaN. + Returns ------- `sunpy.map.Map` diff --git a/sunkit_image/utils/utils.py b/sunkit_image/utils/utils.py index bc506675..b8765062 100644 --- a/sunkit_image/utils/utils.py +++ b/sunkit_image/utils/utils.py @@ -24,6 +24,7 @@ "remove_duplicate", "apply_upsilon", "blackout_pixels_above_radius", + "find_radial_bin_edges" ] @@ -400,3 +401,45 @@ def blackout_pixels_above_radius(smap, radius_limit=1.5 * u.R_sun, fill=np.nan): # Create a new map with the masked data return sunpy.map.Map(masked_data, smap.meta) + + +def find_radial_bin_edges(smap, radial_bin_edges=None): + """ + Calculate radial bin edges for a solar map, either using provided edges or + generating them automatically. + + Parameters + ---------- + smap : `sunpy.map.Map` + A sunpy Map containing the data to be binned. + radial_bin_edges : `astropy.units.Quantity`, optional + Pre-defined bin edges for radial binning. Should be a Quantity array with units + of solar radii (u.R_sun) or pixels. If `None` (the default), bin edges + will be automatically generated based on the map dimensions. + + Returns + ------- + `astropy.units.Quantity` + The final bin edges used for radial binning. + `astropy.units.Quantity` + Array of radial distances for each pixel in the map, matching the input + map dimensions. + """ + # 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