Implemented the Fourier Normalizing Radial Gradient Filter

changelog/17.feature.rst

@@ -0,0 +1 @@
+Added a new feature ```fourier_normalizing_radial_gradient_filter``` to normalize the radial brightness gradient using a fourier approximation.

sunkit_image/offlimb_enhance.py

@@ -6,6 +6,8 @@
 
 import sunpy.map
 
+from sunpy.coordinates import frames
+
 from sunkit_image.utils.utils import (
     bin_edge_summary,
     find_pixel_radii,
@@ -206,10 +208,10 @@ def normalizing_radial_gradient_filter(
     smap,
     radial_bin_edges,
     scale=None,
-    intensity_summary=np.mean,
-    intensity_summary_kwargs=None,
+    intensity_summary=np.nanmean,
+    intensity_summary_kwargs={},
     width_function=np.std,
-    width_function_kwargs=None,
+    width_function_kwargs={},
     application_radius=1 * u.R_sun,
 ):
     """
@@ -232,12 +234,9 @@ def normalizing_radial_gradient_filter(
         For example, in typical Helioprojective Cartesian maps the solar radius is expressed in
         units of arcseconds.
         Defaults to None, which means that the map scale is used.
-    summarize_bin_edges : `str`, optional
-        How to summarize the bin edges.s
-        Defaults to "center".
     intensity_summary : `function`, optional
         A function that returns a summary statistic of the radial intensity.
-        Defaults to `numpy.mean`.
+        Defaults to `numpy.nanmean`.
     intensity_summary_kwargs : None, `dict`
         Keywords applicable to the summary function.
     width_function : `function`
@@ -257,7 +256,8 @@ def normalizing_radial_gradient_filter(
 
     References
     ----------
-    * Morgan, Habbal & Woo, 2006, Sol. Phys., 236, 263. https://link.springer.com/article/10.1007%2Fs11207-006-0113-6
+    * Morgan, Habbal & Woo, 2006, Sol. Phys., 236, 263.
+      https://link.springer.com/article/10.1007%2Fs11207-006-0113-6
     """
 
     # Get the radii for every pixel
@@ -295,3 +295,181 @@ def normalizing_radial_gradient_filter(
         ) / radial_intensity_distribution_summary[i]
 
     return sunpy.map.Map(data, smap.meta)
+
+
+def fourier_normalizing_radial_gradient_filter(
+    smap,
+    radial_bin_edges,
+    order,
+    attenuation_coefficients,
+    ratio_mix=[15, 1],
+    scale=None,
+    intensity_summary=np.nanmean,
+    intensity_summary_kwargs={},
+    width_function=np.std,
+    width_function_kwargs={},
+    application_radius=1 * u.R_sun,
+    number_angular_segments=130,
+):
+    """
+    Implementation of the fourier normalizing radial gradient filter (FNRGF).
+
+    ..note ::
+
+        After applying the filter, current plot settings such as the image normalization
+        may have to be changed in order to obtain a good-looking plot
+
+    Parameters
+    ----------
+    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`
+        Order (Number) of fourier coefficients.
+    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.
+    ratio_mix : `float`
+        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.
+        Defaults to [15,1].
+    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 units
+        of arcseconds.  If None, then the map scale is used.
+    intensity_summary :`function`, optional
+        A function that returns a summary statistic of the radial intensity,
+        Default is numpy.nanmean
+    intensity_summary_kwargs : None, `~dict`
+        Keywords applicable to the summary function.
+    width_function : `function`
+        A function that returns a summary statistic of the distribution of intensity,
+        at a given radius.
+        Defaults to `numpy.std`.
+    width_function_kwargs : `function`
+        Keywords applicable to the width function.
+    application_radius : `astropy.units.Quantity`
+        The FNRGF is applied to emission at radii above the application_radius.
+        Defaults to 1 solar radii.
+    number_angular_segments : `int`
+        Number of angular segments in a circular annulus.
+        Defaults to 130.
+
+    Returns
+    -------
+    new_map : `sunpy.map.Map`
+        A SunPy map that has had the FNRGF applied to it.
+
+    References
+    ----------
+    * Morgan, Habbal & Druckmüllerová, 2011, Astrophysical Journal 737, 88.
+      https://iopscience.iop.org/article/10.1088/0004-637X/737/2/88/pdf.
+    * The implementation is highly inspired by this doctoral thesis.
+      DRUCKMÜLLEROVÁ, H. Application of adaptive filters in processing of solar corona images.
+      https://dspace.vutbr.cz/bitstream/handle/11012/34520/DoctoralThesis.pdf.
+    """
+
+    # Get the radii for every pixel
+    map_r = find_pixel_radii(smap).to(u.R_sun)
+
+    # Get the Helioprojective coordinates of each pixel
+    x, y = np.meshgrid(*[np.arange(v.value) for v in smap.dimensions]) * u.pix
+    coords = smap.pixel_to_world(x, y).transform_to(frames.Helioprojective)
+
+    # Get angles associated with every pixel
+    angles = np.arctan2(coords.Ty.value, coords.Tx.value)
+
+    # Making sure all angles are between (0, 2 * pi)
+    angles = np.where(angles < 0, angles + (2*np.pi), angles)
+
+    # Number of radial bins
+    nbins = radial_bin_edges.shape[1]
+
+    # Storage for the filtered data
+    data = np.zeros_like(smap.data)
+
+    # Iterate over each circular ring
+    for i in range(0, nbins):
+
+        # Finding the pixels which belong to a certain circular ring
+        annulus = np.logical_and(map_r > radial_bin_edges[0, i], map_r < radial_bin_edges[1, i])
+        annulus = np.logical_and(annulus, map_r > application_radius)
+
+        # The angle subtended by each segment
+        segment_angle = 2*np.pi / number_angular_segments
+
+        # Storage of mean and standard deviation of each segment in a circular ring
+        average_segments = np.zeros((1, number_angular_segments))
+        std_dev = np.zeros((1, number_angular_segments))
+
+        # Calculating sin and cos of the angles to be multiplied with the means and standard
+        # deviations to give the fourier approximation
+        cos_matrix = np.cos(np.array([[(2 * np.pi * (j + 1) * (i + 0.5)) / number_angular_segments
+                            for j in range(order)] for i in range(number_angular_segments)]))
+        sin_matrix = np.sin(np.array([[(2 * np.pi * (j + 1) * (i + 0.5)) / number_angular_segments
+                            for j in range(order)] for i in range(number_angular_segments)]))
+
+        # Iterate over each segment in a circular ring
+        for j in range(0, number_angular_segments, 1):
+
+            # Finding all the pixels whose angle values lie in the segment
+            angular_segment = np.logical_and(angles > segment_angle * j,
+                                             angles < segment_angle * (j + 1))
+
+            # Finding the particular segment in the circular ring
+            annulus_segment = np.logical_and(annulus, angular_segment)
+
+            # Finding mean and standard deviation in each segnment. If takes care of the empty
+            # slices.
+            if np.sum([annulus_segment > 0]) == 0:
+                average_segments[0, j] = 0
+                std_dev[0, j] = 0
+            else:
+                average_segments[0, j] = intensity_summary(smap.data[annulus_segment])
+                std_dev[0, j] = width_function(smap.data[annulus_segment])
+
+        # Calculating the fourier coefficients multiplied with the attenuation coefficients
+        fourier_coefficient_a_0 = np.sum(average_segments) * (2 / number_angular_segments)
+        fourier_coefficient_a_0 = fourier_coefficient_a_0 * (attenuation_coefficients[0, 1])
+
+        fourier_coefficients_a_k = np.matmul(average_segments, cos_matrix) * (2 / number_angular_segments)
+        fourier_coefficients_a_k = fourier_coefficients_a_k * attenuation_coefficients[0][1:]
+
+        fourier_coefficients_b_k = np.matmul(average_segments, sin_matrix) * (2 / number_angular_segments)
+        fourier_coefficients_b_k = fourier_coefficients_b_k * attenuation_coefficients[0][1:]
+
+        fourier_coefficient_c_0 = np.sum(std_dev) * (2 / number_angular_segments)
+        fourier_coefficient_c_0 = fourier_coefficient_c_0 * attenuation_coefficients[1, 1]
+
+        fourier_coefficients_c_k = np.matmul(std_dev, cos_matrix) * (2 / number_angular_segments)
+        fourier_coefficients_c_k = fourier_coefficients_c_k * attenuation_coefficients[1][1:]
+
+        fourier_coefficients_d_k = np.matmul(std_dev, sin_matrix) * (2 / number_angular_segments)
+        fourier_coefficients_d_k = fourier_coefficients_d_k * attenuation_coefficients[1][1:]
+
+        # To calculate the multiples of angles of each pixel for finding the fourier approximation
+        # at that point. See equations 6.8 and 6.9 of the doctoral thesis.
+        K_matrix = np.ones((order, np.sum(annulus > 0))) * np.array(range(1, order+1)).T.reshape(order, 1)
+        phi_matrix = angles[annulus].reshape((1, angles[annulus].shape[0]))
+        angles_of_pixel = K_matrix * phi_matrix
+
+        # Get the approxiamted value of mean
+        mean_approximated = np.matmul(fourier_coefficients_a_k, np.cos(angles_of_pixel))
+        mean_approximated = mean_approximated + np.matmul(fourier_coefficients_b_k, np.sin(angles_of_pixel))
+        mean_approximated = mean_approximated + fourier_coefficient_a_0 / 2
+
+        # Get the approxiamted value of standard deviation
+        std_approximated = np.matmul(fourier_coefficients_c_k, np.cos(angles_of_pixel))
+        std_approximated = std_approximated + np.matmul(fourier_coefficients_d_k, np.sin(angles_of_pixel))
+        std_approximated = std_approximated + fourier_coefficient_c_0 / 2
+
+        # Normailize the data
+        data[annulus] = np.ravel((smap.data[annulus] - mean_approximated) / std_approximated)
+
+        # Linear combination of original image and the filtered data.
+        data[annulus] = ratio_mix[0] * smap.data[annulus] + ratio_mix[1] * data[annulus]
+
+    return sunpy.map.Map(data, smap.meta)

sunkit_image/tests/test_offlimb_enhance.py

@@ -1,12 +1,55 @@
-import os
-
 import numpy as np
 import pytest
 import astropy.units as u
-from numpy.testing import assert_allclose
 
 import sunpy
 import sunpy.map
 import sunpy.data.test
+import sunkit_image.utils.utils as utils
+import sunkit_image.offlimb_enhance as off
+
+
+test_map_data = np.ones((5, 5))
+header = {}
+test_map = sunpy.map.Map((test_map_data, header))
+radial_bin_edges = utils.equally_spaced_bins()
+radial_bin_edges = radial_bin_edges*u.R_sun
+
+
+def set_attenuation_coefficients(order):
+    attenuation_coefficients = np.zeros((2, order + 1))
+    attenuation_coefficients[0][:] = np.linspace(1, 0, order + 1)
+    attenuation_coefficients[1][:] = np.linspace(1, 0, order + 1)
+    return attenuation_coefficients
+
+
+def test_normalizing_radial_gradient_filter():
+
+    with pytest.warns(Warning, match="Missing metadata for solar radius: assuming photospheric limb as seen from Earth"):
+
+        result = np.zeros_like(test_map_data)
+        expect = off.normalizing_radial_gradient_filter(test_map, radial_bin_edges)
+
+        assert np.allclose(expect.data.shape, test_map.data.shape)
+        assert np.allclose(expect.data, result)
+
+
+def test_fourier_normalizing_radial_gradient_filter():
+
+    with pytest.warns(Warning, match="Missing metadata for solar radius: assuming photospheric limb as seen from Earth"):
+
+        result = np.zeros_like(test_map_data)
+
+        order = 20
+        attenuation_coefficients = set_attenuation_coefficients(order)
+        expect = off.fourier_normalizing_radial_gradient_filter(test_map, radial_bin_edges, order,
+                                                                attenuation_coefficients)
+        assert np.allclose(expect.data.shape, test_map.data.shape)
+        assert np.allclose(expect.data, result)
 
-# TODO: Write this
+        order = 33
+        attenuation_coefficients = set_attenuation_coefficients(order)
+        expect = off.fourier_normalizing_radial_gradient_filter(test_map, radial_bin_edges, order,
+                                                                attenuation_coefficients)
+        assert np.allclose(expect.data.shape, test_map.data.shape)
+        assert np.allclose(expect.data, result)