Skip to content

Commit

Permalink
Update code
Browse files Browse the repository at this point in the history
  • Loading branch information
nabobalis committed Oct 29, 2024
1 parent 9dd000d commit cd8f0eb
Show file tree
Hide file tree
Showing 11 changed files with 215 additions and 168 deletions.
1 change: 1 addition & 0 deletions changelog/205.doc.rst
Original file line number Diff line number Diff line change
@@ -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.
1 change: 1 addition & 0 deletions changelog/231.breaking.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
All functions under (`sunkit_image.radial`) now have a settable fill value which now defaults to NaN instead of 0.
1 change: 1 addition & 0 deletions changelog/231.bugfix.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Improved the Radial Histogram Equalizing Filter (RHEF) (`sunkit_image.radial.rhef`) to work with images outside of SDO/AIA.
6 changes: 3 additions & 3 deletions examples/radial_gradient_filters.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@
radial_bin_edges = equally_spaced_bins(1, 2, aia_map.data.shape[0] // 4)
radial_bin_edges *= u.R_sun

base_nrgf = radial.nrgf(aia_map, radial_bin_edges, application_radius=1 * u.R_sun, progress=False)
base_nrgf = radial.nrgf(aia_map, radial_bin_edges=radial_bin_edges, application_radius=1 * u.R_sun)

###########################################################################
# We will need to work out a few parameters for the FNRGF.
Expand All @@ -50,12 +50,12 @@
order = 20
attenuation_coefficients = radial.set_attenuation_coefficients(order)

base_fnrgf = radial.fnrgf(aia_map, radial_bin_edges, order, attenuation_coefficients, application_radius=1 * u.R_sun, progress=False)
base_fnrgf = radial.fnrgf(aia_map, radial_bin_edges=radial_bin_edges, order=order, attenuation_coefficients=attenuation_coefficients, application_radius=1 * u.R_sun)

###########################################################################
# Now we will also use the final filter, RHEF.

base_rhef = radial.rhef(aia_map, radial_bin_edges, application_radius=1 * u.R_sun, progress=False)
base_rhef = radial.rhef(aia_map, radial_bin_edges=radial_bin_edges, application_radius=1 * u.R_sun)

###########################################################################
# Finally we will plot the filtered maps with the original to demonstrate the effect of each.
Expand Down
8 changes: 7 additions & 1 deletion pytest.ini
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,12 @@ mpl-use-full-test-name = True
addopts =
--doctest-rst
-p no:unraisableexception
-p no:threadexception
-p no:theadexception
-m "not mpl_image_compare"
--dist no
--arraydiff
--doctest-ignore-import-errors
--doctest-continue-on-failure
markers =
remote_data: marks this test function as needing remote data.
online: marks this test function as needing online connectivity.
Expand Down Expand Up @@ -59,3 +64,4 @@ filterwarnings =
# astropy/utils/iers/iers.py:1115: in auto_open good_enough = cls._today() + TimeDelta(offset, format="jd")
ignore:datetime.datetime.utc:DeprecationWarning
ignore:leap-second file is expired:astropy.utils.iers.iers.IERSStaleWarning
ignore:leap-second auto-update failed:astropy.utils.exceptions.AstropyWarning
18 changes: 16 additions & 2 deletions sunkit_image/conftest.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
import os
import logging
import tempfile
import warnings
import importlib.util

import numpy as np
Expand All @@ -9,6 +10,7 @@

import astropy
import astropy.config.paths
from astropy.io.fits.verify import VerifyWarning
from astropy.utils.data import get_pkg_data_filename

import sunpy.data.sample
Expand Down Expand Up @@ -170,15 +172,27 @@ def granule_minimap3():
)
return sunpy.map.GenericMap(arr, header)


@pytest.fixture(params=["array", "map"])
def aia_171(request):
smap = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE)
# VerifyWarning: Invalid 'BLANK' keyword in header.
# The 'BLANK' keyword is only applicable to integer data, and will be ignored in this HDU.
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=VerifyWarning)
smap = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE)
if request.param == "map":
return smap
return smap.data if request.param == "array" else None


@pytest.fixture()
def aia_171_map():
# VerifyWarning: Invalid 'BLANK' keyword in header.
# The 'BLANK' keyword is only applicable to integer data, and will be ignored in this HDU.
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=VerifyWarning)
return sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE)


@pytest.fixture()
def hmi_map():
hmi_file = get_test_filepath("hmi_continuum_test_lowres_data.fits")
Expand Down
160 changes: 75 additions & 85 deletions sunkit_image/radial.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
)

Expand Down Expand Up @@ -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

Check warning on line 106 in sunkit_image/radial.py

View check run for this annotation

Codecov / codecov/patch

sunkit_image/radial.py#L105-L106

Added lines #L105 - L106 were not covered by tests

return stats.rankdata(arr, method="average") / len(arr)

Check warning on line 108 in sunkit_image/radial.py

View check run for this annotation

Codecov / codecov/patch

sunkit_image/radial.py#L108

Added line #L108 was not covered by tests

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))

Check warning on line 114 in sunkit_image/radial.py

View check run for this annotation

Codecov / codecov/patch

sunkit_image/radial.py#L110-L114

Added lines #L110 - L114 were not covered by tests

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))

Check warning on line 119 in sunkit_image/radial.py

View check run for this annotation

Codecov / codecov/patch

sunkit_image/radial.py#L116-L119

Added lines #L116 - L119 were not covered by tests

# 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

Check warning on line 127 in sunkit_image/radial.py

View check run for this annotation

Codecov / codecov/patch

sunkit_image/radial.py#L122-L127

Added lines #L122 - L127 were not covered by tests
else:
msg = f"{method} is invalid. Allowed values are 'inplace', 'numpy', 'scipy'"

Check warning on line 129 in sunkit_image/radial.py

View check run for this annotation

Codecov / codecov/patch

sunkit_image/radial.py#L129

Added line #L129 was not covered by tests
raise NotImplementedError(msg)
return ranking_func

Check warning on line 131 in sunkit_image/radial.py

View check run for this annotation

Codecov / codecov/patch

sunkit_image/radial.py#L131

Added line #L131 was not covered by tests

def intensity_enhance(
smap,
*,
radial_bin_edges=None,
scale=None,
summarize_bin_edges="center",
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -209,22 +238,22 @@ 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


def nrgf(
smap,
*,
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,
progress=False,
fill=np.nan,
):
"""
Expand All @@ -247,9 +276,10 @@ def nrgf(
----------
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 : None or `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
Expand All @@ -269,11 +299,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
-------
Expand Down Expand Up @@ -394,15 +425,16 @@ def set_attenuation_coefficients(order, range_mean=None, range_std=None, cutoff=

def fnrgf(
smap,
radial_bin_edges,
order,
*,
attenuation_coefficients,
radial_bin_edges=None,
order=3,
ratio_mix=None,
intensity_summary=np.nanmean,
width_function=np.std,
application_radius=1 * u.R_sun,
number_angular_segments=130,
progress=True,
progress=False,
fill=np.nan,
):
"""
Expand All @@ -427,10 +459,13 @@ def fnrgf(
----------
smap : `sunpy.map.Map`
A SunPy map.
radial_bin_edges : `astropy.units.Quantity`
A two-dimensional array of bin edges of size ``[2, nbins]`` where ``nbins`` is the number of bins.
order : `int`
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.
order : `int`, optional
Order (number) of fourier coefficients and it can not be lower than 1.
Defaults to 3.
attenuation_coefficients : `float`
A two dimensional array of shape ``[2, order + 1]``. The first row contain attenuation
coefficients for mean calculations. The second row contains attenuation coefficients
Expand All @@ -451,11 +486,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
-------
Expand Down Expand Up @@ -595,57 +631,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,
Expand All @@ -655,7 +640,7 @@ def rhef(
upsilon=0.35,
method="numpy",
vignette=None,
progress=True,
progress=False,
fill=np.nan,
):
"""
Expand All @@ -677,21 +662,27 @@ def rhef(
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.
These define the radial segments where filtering is applied.
If None, radial bins will be generated automatically.
application_radius : `astropy.units.Quantity`, optional
The radius above which to apply the RHEF. Only regions with radii above this value will be filtered.
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.
Radius beyond which pixels will be set to NaN.
Must be in units that are compatible with "R_sun" as the value will be transformed.
Defaults to `None`.
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`
Expand Down Expand Up @@ -734,7 +725,6 @@ def rhef(
# Adjust plot settings to remove extra normalization
# 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
Expand Down
Loading

0 comments on commit cd8f0eb

Please sign in to comment.