diff --git a/sunkit_image/radial.py b/sunkit_image/radial.py index 9a460ae5..a7b0623e 100644 --- a/sunkit_image/radial.py +++ b/sunkit_image/radial.py @@ -15,17 +15,16 @@ apply_upsilon, bin_edge_summary, blackout_pixels_above_radius, - equally_spaced_bins, - find_pixel_radii, + find_radial_bin_edges, get_radial_intensity_summary, ) -__all__ = ["fnrgf", "intensity_enhance", "set_attenuation_coefficients", "nrgf", "rhef"] +__all__ = ["fnrgf", "intensity_enhance", "nrgf", "rhef"] def _fit_polynomial_to_log_radial_intensity(radii, intensity, degree): """ - Fits a polynomial of a given degree to the log of the radial intensity. + Fits a polynomial of a given degree to the log of the radial intensity. Parameters ---------- @@ -102,9 +101,41 @@ def _normalize_fit_radial_intensity(radii, polynomial, normalization_radius): ) +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, + *, + radial_bin_edges=None, scale=None, summarize_bin_edges="center", summary=np.mean, @@ -135,9 +166,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 @@ -169,8 +201,8 @@ def intensity_enhance( `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,6 +239,7 @@ def intensity_enhance( 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 @@ -214,14 +247,16 @@ def intensity_enhance( 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, + progress=False, + fill=np.nan, ): """ Implementation of the normalizing radial gradient filter (NRGF). @@ -243,9 +278,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 @@ -266,8 +302,11 @@ def nrgf( 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. + 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 ------- @@ -285,15 +324,9 @@ def nrgf( 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 +347,7 @@ def nrgf( ) # 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): @@ -329,80 +362,21 @@ def nrgf( return new_map -def set_attenuation_coefficients(order, range_mean=None, range_std=None, cutoff=0): - """ - This is a helper function to Fourier Normalizing Radial Gradient Filter - (`sunkit_image.radial.fnrgf`). - - This function sets the attenuation coefficients in the one of the following two manners: - - If ``cutoff`` is ``0``, then it will set the attenuation coefficients as linearly decreasing between - the range ``range_mean`` for the attenuation coefficients for mean approximation and ``range_std`` for - the attenuation coefficients for standard deviation approximation. - - If ``cutoff`` is not ``0``, then it will set the last ``cutoff`` number of coefficients equal to zero - while all the others the will be set as linearly decreasing as described above. - - .. note:: - - This function only describes some of the ways in which attenuation coefficients can be calculated. - The optimal coefficients depends on the size and quality of image. There is no generalized formula - for choosing them and its up to the user to choose a optimum value. - - .. note:: - - The returned maps have their ``plot_settings`` changed to remove the extra normalization step. - - Parameters - ---------- - order : `int` - The order of the Fourier approximation. - range_mean : `list`, optional - A list of length of ``2`` which contains the highest and lowest values between which the coefficients for - mean approximation be calculated in a linearly decreasing manner. - range_std : `list`, optional - A list of length of ``2`` which contains the highest and lowest values between which the coefficients for - standard deviation approximation be calculated in a linearly decreasing manner. - cutoff : `int`, optional - The numbers of coefficients from the last that should be set to ``zero``. - - Returns - ------- - `numpy.ndarray` - A numpy array of shape ``[2, order + 1]`` containing the attenuation coefficients for the Fourier - coffiecients. The first row describes the attenustion coefficients for the Fourier coefficients of - the mean approximation. The second row contains the attenuation coefficients for the Fourier coefficients - of the standard deviation approximation. - """ - if range_std is None: - range_std = [1.0, 0.0] - if range_mean is None: - range_mean = [1.0, 0.0] - attenuation_coefficients = np.zeros((2, order + 1)) - attenuation_coefficients[0, :] = np.linspace(range_mean[0], range_mean[1], order + 1) - attenuation_coefficients[1, :] = np.linspace(range_std[0], range_std[1], order + 1) - - if cutoff > (order + 1): - msg = "Cutoff cannot be greater than order + 1." - raise ValueError(msg) - - if cutoff != 0: - attenuation_coefficients[:, (-1 * cutoff) :] = 0 - - return attenuation_coefficients - - def fnrgf( smap, - radial_bin_edges, - order, - attenuation_coefficients, + *, + radial_bin_edges=None, + order=3, + range_mean=None, + range_std=None, + cutoff=0, 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, ): """ Implementation of the fourier normalizing radial gradient filter (FNRGF). @@ -426,14 +400,21 @@ 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. - 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 - for standard deviation calculation. + Defaults to 3. + range_mean : `list`, optional + A list of length of ``2`` which contains the highest and lowest values between which the coefficients for + mean approximation be calculated in a linearly decreasing manner. + range_std : `list`, optional + A list of length of ``2`` which contains the highest and lowest values between which the coefficients for + standard deviation approximation be calculated in a linearly decreasing manner. + cutoff : `int`, optional + The numbers of coefficients from the last that should be set to ``zero``. ratio_mix : `float`, optional A one dimensional array of shape ``[2, 1]`` with values equal to ``[K1, K2]``. The ratio in which the original image and filtered image are mixed. @@ -451,8 +432,11 @@ def fnrgf( Number of angular segments in a circular annulus. Defaults to 130. progress : `bool`, optional - Display a progressbar on the main loop. - Defaults to True. + 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 ------- @@ -474,16 +458,8 @@ def fnrgf( 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 +475,23 @@ def fnrgf( nbins = radial_bin_edges.shape[1] # Storage for the filtered data - data = np.zeros_like(smap.data) + data = np.ones_like(smap.data) * fill + + # Set attenuation coefficients + if range_std is None: + range_std = [1.0, 0.0] + if range_mean is None: + range_mean = [1.0, 0.0] + attenuation_coefficients = np.zeros((2, order + 1)) + attenuation_coefficients[0, :] = np.linspace(range_mean[0], range_mean[1], order + 1) + attenuation_coefficients[1, :] = np.linspace(range_std[0], range_std[1], order + 1) + + if cutoff > (order + 1): + msg = "Cutoff cannot be greater than order + 1." + raise ValueError(msg) + + if cutoff != 0: + attenuation_coefficients[:, (-1 * cutoff) :] = 0 # Iterate over each circular ring for i in tqdm(range(nbins), desc="FNRGF: ", disable=not progress): @@ -600,47 +592,17 @@ 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 - - @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, - progress=True, + vignette=None, + progress=False, + fill=np.nan, ): """ Implementation of the Radial Histogram Equalizing Filter (RHEF). @@ -658,33 +620,34 @@ def rhef( 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 : ``{"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. + 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` - 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 +658,35 @@ def rhef( https://www.proquest.com/docview/2759080511 """ - # Get the radii for every pixel - map_r = find_pixel_radii(smap).to(u.R_sun) - - 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 - - # 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], - ) + radial_bin_edges, map_r = find_radial_bin_edges(smap, radial_bin_edges) - # Allocate storage for the filtered data - data = np.zeros_like(smap.data) - meta = smap.meta + data = np.ones_like(smap.data) * fill - 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)) - # This must be done whenever one is adjusting the overall statistical distribution of values + # 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 return new_map diff --git a/sunkit_image/tests/test_radial.py b/sunkit_image/tests/test_radial.py index 9f38fd2c..d8b4c95d 100644 --- a/sunkit_image/tests/test_radial.py +++ b/sunkit_image/tests/test_radial.py @@ -49,13 +49,13 @@ def radial_bin_edges(): def test_nrgf(map_test1, map_test2, radial_bin_edges): result = np.zeros_like(map_test1.data) - expect = rad.nrgf(map_test1, radial_bin_edges, application_radius=0.001 * u.R_sun) + expect = rad.nrgf(map_test1, radial_bin_edges=radial_bin_edges, application_radius=0.001 * u.R_sun, fill=0) assert np.allclose(expect.data.shape, map_test1.data.shape) assert np.allclose(expect.data, result) # Hand calculated - result = [ + result1 = [ [0.0, 1.0, 0.0, -1.0, 0.0], [1.0, 0.0, 0.0, 0.0, -1.0], [0.0, 0.0, 0.0, 0.0, 0.0], @@ -63,16 +63,14 @@ def test_nrgf(map_test1, map_test2, radial_bin_edges): [0.0, 1.0, 0.0, -1.0, 0.0], ] - expect = rad.nrgf(map_test2, radial_bin_edges, application_radius=0.001 * u.R_sun) + expect1 = rad.nrgf(map_test2, radial_bin_edges=radial_bin_edges, application_radius=0.001 * u.R_sun, fill=0) - assert np.allclose(expect.data.shape, map_test2.data.shape) - assert np.allclose(expect.data, result) + assert np.allclose(expect1.data.shape, map_test2.data.shape) + assert np.allclose(expect1.data, result1) def test_fnrgf(map_test1, map_test2, radial_bin_edges): order = 1 - - # Hand calculated result = [ [-0.0, 96.0, 128.0, 96.0, -0.0], [96.0, 224.0, 288.0, 224.0, 96.0], @@ -80,44 +78,44 @@ def test_fnrgf(map_test1, map_test2, radial_bin_edges): [96.0, 224.0, 288.0, 224.0, 96.0], [-0.0, 96.0, 128.0, 96.0, -0.0], ] - attenuation_coefficients = rad.set_attenuation_coefficients(order) expect = rad.fnrgf( map_test1, - radial_bin_edges, - order, - attenuation_coefficients, + radial_bin_edges=radial_bin_edges, + order=order, + range_mean=[1.0, 0.0], + range_std=[1.0, 0.0], + cutoff=0, application_radius=0.001 * u.R_sun, number_angular_segments=4, + fill=0, ) assert np.allclose(expect.data.shape, map_test1.data.shape) assert np.allclose(expect.data, result) - # Hand calculated - result = [ + result1 = [ [-0.0, 128.0, 128.0, 96.0, -0.0], [128.0, 224.0, 288.0, 224.0, 96.0], [128.0, 288.0, 0.0, 288.0, 128.0], [128.0, 224.0, 288.0, 224.0, 96.0], [-0.0, 128.0, 128.0, 96.0, -0.0], ] - expect = rad.fnrgf( + expect1 = rad.fnrgf( map_test2, - radial_bin_edges, - order, - attenuation_coefficients, + radial_bin_edges=radial_bin_edges, + order=order, + range_mean=[1.0, 0.0], + range_std=[1.0, 0.0], + cutoff=0, application_radius=0.001 * u.R_sun, number_angular_segments=4, + fill=0, ) - assert np.allclose(expect.data.shape, map_test2.data.shape) - assert np.allclose(expect.data, result) + assert np.allclose(expect1.data.shape, map_test2.data.shape) + assert np.allclose(expect1.data, result1) - # The below tests are dummy tests. - # These values were not verified by hand rather they were - # generated using the code itself. order = 5 - - result = [ + result2 = [ [-0.0, 90.52799999982116, 126.73137084989847, 90.52799999984676, -0.0], [90.52800000024544, 207.2, 285.14558441227155, 207.2, 90.5280000001332], [126.73137084983244, 285.1455844119744, 0.0, 280.05441558770406, 124.4686291500961], @@ -125,20 +123,22 @@ def test_fnrgf(map_test1, map_test2, radial_bin_edges): [0.0, 90.52799999986772, 124.46862915010152, 90.52799999989331, -0.0], ] - attenuation_coefficients = rad.set_attenuation_coefficients(order) - expect = rad.fnrgf( + expect2 = rad.fnrgf( map_test1, - radial_bin_edges, - order, - attenuation_coefficients, + radial_bin_edges=radial_bin_edges, + order=order, + range_mean=[1.0, 0.0], + range_std=[1.0, 0.0], + cutoff=0, application_radius=0.001 * u.R_sun, number_angular_segments=4, + fill=0, ) - assert np.allclose(expect.data.shape, map_test1.data.shape) - assert np.allclose(expect.data, result) + assert np.allclose(expect2.data.shape, map_test1.data.shape) + assert np.allclose(expect2.data, result2) - result = [ + result3 = [ [-0.0, 120.55347470594926, 126.73137084989847, 90.67852529365966, -0.0], [120.70526403418884, 207.2, 285.14558441227155, 207.2, 90.52673596626707], [126.73137084983244, 285.1455844119744, 0.0, 280.05441558770406, 124.4686291500961], @@ -146,73 +146,64 @@ def test_fnrgf(map_test1, map_test2, radial_bin_edges): [0.0, 120.55347470601022, 124.46862915010152, 90.67852529370734, -0.0], ] - attenuation_coefficients = rad.set_attenuation_coefficients(order) - expect = rad.fnrgf( + expect3 = rad.fnrgf( map_test2, - radial_bin_edges, - order, - attenuation_coefficients, + radial_bin_edges=radial_bin_edges, + order=order, + range_mean=[1.0, 0.0], + range_std=[1.0, 0.0], + cutoff=0, application_radius=0.001 * u.R_sun, number_angular_segments=4, + fill=0, ) - assert np.allclose(expect.data.shape, map_test2.data.shape) - assert np.allclose(expect.data, result) + assert np.allclose(expect3.data.shape, map_test2.data.shape) + assert np.allclose(expect3.data, result3) - order = 0 +def test_fnrgf_errors(map_test1): with pytest.raises(ValueError, match="Minimum value of order is 1"): rad.fnrgf( - map_test2, - radial_bin_edges, - order, - attenuation_coefficients, - application_radius=0.001 * u.R_sun, - number_angular_segments=4, + map_test1, + order=0, + range_mean=[1.0, 0.0], + range_std=[1.0, 0.0], + cutoff=0, ) - -@pytest.fixture() -def smap(): - return sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE) - - @figure_test @pytest.mark.remote_data() -def test_fig_nrgf(smap): +def test_fig_nrgf(aia_171_map): radial_bin_edges = utils.equally_spaced_bins() radial_bin_edges *= u.R_sun - out = rad.nrgf(smap, radial_bin_edges) - + out = rad.nrgf(aia_171_map, radial_bin_edges=radial_bin_edges) out.plot() @figure_test @pytest.mark.remote_data() -def test_fig_fnrgf(smap): +def test_fig_fnrgf(aia_171_map): radial_bin_edges = utils.equally_spaced_bins() radial_bin_edges *= u.R_sun - order = 20 - attenuation_coefficients = rad.set_attenuation_coefficients(order) - out = rad.fnrgf(smap, radial_bin_edges, order, attenuation_coefficients) - + out = rad.fnrgf(aia_171_map, radial_bin_edges=radial_bin_edges, order=order, range_mean=[1.0, 0.0], range_std=[1.0, 0.0], cutoff=0) out.plot() @figure_test @pytest.mark.remote_data() -def test_fig_rhef(smap): - radial_bin_edges = utils.equally_spaced_bins(0, 2, smap.data.shape[1]) +def test_fig_rhef(aia_171_map): + radial_bin_edges = utils.equally_spaced_bins(0, 2, aia_171_map.data.shape[1]) radial_bin_edges *= u.R_sun - out = rad.rhef(smap, radial_bin_edges=radial_bin_edges, upsilon=0.35, method="scipy") + out = rad.rhef(aia_171_map, radial_bin_edges=radial_bin_edges, upsilon=None, method="scipy") out.plot() @figure_test @pytest.mark.remote_data() -def test_multifig_rhef(smap): - radial_bin_edges = utils.equally_spaced_bins(0, 2, smap.data.shape[1]) +def test_multifig_rhef(aia_171_map): + radial_bin_edges = utils.equally_spaced_bins(0, 2, aia_171_map.data.shape[1]) radial_bin_edges *= u.R_sun # Define the list of upsilon pairs where the first number affects dark components and the second number affects bright ones @@ -225,9 +216,9 @@ def test_multifig_rhef(smap): ] # Crop the figures to see better detail - top_right = SkyCoord(1200 * u.arcsec, 0 * u.arcsec, frame=smap.coordinate_frame) - bottom_left = SkyCoord(0 * u.arcsec, -1200 * u.arcsec, frame=smap.coordinate_frame) - aia_map_cropped = smap.submap(bottom_left, top_right=top_right) + top_right = SkyCoord(1200 * u.arcsec, 0 * u.arcsec, frame=aia_171_map.coordinate_frame) + bottom_left = SkyCoord(0 * u.arcsec, -1200 * u.arcsec, frame=aia_171_map.coordinate_frame) + aia_map_cropped = aia_171_map.submap(bottom_left, top_right=top_right) fig, axes = plt.subplots( 2, 3, figsize=(15, 10), sharex="all", sharey="all", subplot_kw={"projection": aia_map_cropped} ) @@ -238,7 +229,7 @@ def test_multifig_rhef(smap): # Loop through the upsilon_list and plot each filtered map for i, upsilon in enumerate(upsilon_list): - out_map = rad.rhef(smap, upsilon=upsilon, method="scipy") + out_map = rad.rhef(aia_171_map, upsilon=upsilon, method="scipy") out_map_crop = out_map.submap(bottom_left, top_right=top_right) out_map_crop.plot(axes=axes[i + 1]) axes[i + 1].set_title(f"Upsilon = {upsilon}") @@ -248,30 +239,6 @@ def test_multifig_rhef(smap): return fig -def test_set_attenuation_coefficients(): - order = 1 - # Hand calculated - expect1 = [[1, 0.0], [1, 0.0]] - - result1 = rad.set_attenuation_coefficients(order) - assert np.allclose(expect1, result1) - - order = 3 - # Hand calculated - expect2 = [[1.0, 0.66666667, 0.33333333, 0.0], [1.0, 0.66666667, 0.33333333, 0.0]] - - result2 = rad.set_attenuation_coefficients(order) - assert np.allclose(expect2, result2) - - expect3 = [[1.0, 0.66666667, 0.0, 0.0], [1.0, 0.66666667, 0.0, 0.0]] - - result3 = rad.set_attenuation_coefficients(order, cutoff=2) - assert np.allclose(expect3, result3) - - with pytest.raises(ValueError, match="Cutoff cannot be greater than order \\+ 1"): - rad.set_attenuation_coefficients(order, cutoff=5) - - def test_fit_polynomial_to_log_radial_intensity(): radii = (0.001, 0.002) * u.R_sun intensity = np.asarray([1, 2]) @@ -330,10 +297,15 @@ def test_intensity_enhance(map_test1): enhancement = 1 / rad._normalize_fit_radial_intensity(map_r, polynomial, normalization_radius) # NOQA: SLF001 enhancement[map_r < normalization_radius] = 1 - with pytest.raises(ValueError, match="The fit range must be strictly increasing."): - rad.intensity_enhance(smap=map_test1, radial_bin_edges=radial_bin_edges, scale=scale, fit_range=fit_range[::-1]) - assert np.allclose( enhancement * map_test1.data, - rad.intensity_enhance(smap=map_test1, radial_bin_edges=radial_bin_edges, scale=scale).data, + rad.intensity_enhance(map_test1, radial_bin_edges=radial_bin_edges, scale=scale).data, ) + + +@skip_windows +def test_intensity_enhance_errors(map_test1): + fit_range = [1, 1.5] * u.R_sun + scale = 1 * map_test1.rsun_obs + with pytest.raises(ValueError, match="The fit range must be strictly increasing."): + rad.intensity_enhance(map_test1, scale=scale, fit_range=fit_range[::-1])