Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added an example for advanced usage of WOW #238

Merged
merged 4 commits into from
Oct 29, 2024
Merged
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
96 changes: 96 additions & 0 deletions examples/advanced_wow.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,96 @@
"""
====================================================
Advanced Usage of Wavelets Optimized Whitening (WOW)
====================================================

This example demonstrates different options of the Wavelets Optimized Whitening applied to a `sunpy.map.Map`
using `sunkit_image.enhance.wow`.
"""

import matplotlib.pyplot as plt

import sunpy.data.sample
import sunpy.map

import sunkit_image.enhance as enhance

from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.visualization import ImageNormalize, LinearStretch, PowerStretch, AsymmetricPercentileInterval

###########################################################################
# `sunpy` provides a range of sample data with a number of suitable images.
# Here will just use AIA 193
nabobalis marked this conversation as resolved.
Show resolved Hide resolved

input_map = sunpy.map.Map(sunpy.data.sample.AIA_193_JUN2012)

###########################################################################
# We crop the southwestern quadrant
nabobalis marked this conversation as resolved.
Show resolved Hide resolved

top_right = SkyCoord(1200 * u.arcsec, 0 * u.arcsec, frame=input_map.coordinate_frame)
bottom_left = SkyCoord(0 * u.arcsec, -1200 * u.arcsec, frame=input_map.coordinate_frame)
input_map = input_map.submap(bottom_left, top_right=top_right)

###########################################################################
# We now applying different options of the Wavelets Optimized Whitening alogrithm.
# The `sunkit_image.enhance.wow` function takes either a `sunpy.map.Map` or a `numpy.ndarray` as a input.
# First, we call WOW with now arguments, which returns the default WOW enhancement.
nabobalis marked this conversation as resolved.
Show resolved Hide resolved

wow_map = enhance.wow(input_map)

###########################################################################
# The we denoise the output using soft thresholding in the three first wavelet scales with sigma = 5, 2, 1
nabobalis marked this conversation as resolved.
Show resolved Hide resolved

denoise_coefficients = [5, 2, 1]
wow_map_denoised = enhance.wow(input_map,
denoise_coefficients=denoise_coefficients)
nabobalis marked this conversation as resolved.
Show resolved Hide resolved

###########################################################################
# We then run the edge-aware (bilateral) flavor of the algorithm. Its prevents ringing around sharp edges (e.g. the
# solar limb, of very bright features.

wow_map_bilateral = enhance.wow(input_map,
bilateral=1)
nabobalis marked this conversation as resolved.
Show resolved Hide resolved

###########################################################################
# This will call the edge-aware algorithm with denoising.

wow_map_bilateral_denoised = enhance.wow(input_map,
bilateral=1,
denoise_coefficients=denoise_coefficients)
nabobalis marked this conversation as resolved.
Show resolved Hide resolved

###########################################################################
# Finally, we merge the denoised edge-aware enhanced image with the gamma-stretched input, with weight h.

gamma = 4
wow_map_bilateral_denoised_merged = enhance.wow(input_map,
bilateral=1,
denoise_coefficients=denoise_coefficients,
gamma=gamma,
h=0.99)

nabobalis marked this conversation as resolved.
Show resolved Hide resolved
###########################################################################
# Now we will plot the final result and compare that to the original image.
nabobalis marked this conversation as resolved.
Show resolved Hide resolved

fig = plt.figure(figsize=(8, 12))

nabobalis marked this conversation as resolved.
Show resolved Hide resolved
variations = {
'Input | gamma = {gamma} stretch': {'map': input_map, 'stretch': PowerStretch(1 / gamma)},
'WOW | linear stretch': {'map': wow_map, 'stretch': LinearStretch()},
'denoised WOW': {'map': wow_map_denoised, 'stretch': LinearStretch()},
'Edge-aware WOW': {'map': wow_map_bilateral, 'stretch': LinearStretch()},
'Edge-aware & denoised WOW': {'map': wow_map_bilateral_denoised, 'stretch': LinearStretch()},
'Merged with input': {'map': wow_map_bilateral_denoised_merged, 'stretch': LinearStretch()}
}

nabobalis marked this conversation as resolved.
Show resolved Hide resolved
interval = AsymmetricPercentileInterval(1, 99.9)

nabobalis marked this conversation as resolved.
Show resolved Hide resolved
for i, (title, image) in enumerate(variations.items()):
ax = fig.add_subplot(3, 2, i + 1, projection=image['map'])
image['map'].plot(norm=ImageNormalize(image['map'].data, interval=interval, stretch=image['stretch']))
ax.set_title(title)
ax.axis('off')

fig.tight_layout()

plt.show()