Skip to content

Commit

Permalink
updated examples
Browse files Browse the repository at this point in the history
  • Loading branch information
sfarrens committed Jun 23, 2022
1 parent 6167df3 commit 5f75837
Show file tree
Hide file tree
Showing 2 changed files with 94 additions and 79 deletions.
85 changes: 45 additions & 40 deletions examples/plot_galaxy_deconvolution.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,88 +2,93 @@
Galaxy Image Deconvolution
==========================
Credit: S. Farrens
.. codeauthor:: Samuel Farrens <[email protected]>
In this tutorial we will deconvolve the PSF effects from an example galaxy
image.
image using the PySAP-Astro plug-in.
Import Dependencies
Import dependencies
-------------------
Import functions from PySAP and ModOpt.
Import functions from PySAP and
`ModOpt <https://cea-cosmic.github.io/ModOpt/>`_.
"""

import numpy as np
import matplotlib.pyplot as plt
from pysap import Image
from pysap.data import get_sample_data
from astro.deconvolution.deconvolve import sparse_deconv_condatvu
from modopt.signal.noise import add_noise
from modopt.math.convolve import convolve

#############################################################################
# %%
# Load astro data
# ---------------
#
# Load the example images
# First, we load some example images from the PySAP sample data sets.

galaxy = get_sample_data('astro-galaxy')
psf = get_sample_data('astro-psf')

#############################################################################
# Show the clean galaxy image
# %%
# Then we can show the clean galaxy image that we will attempt to recover.

# galaxy.show()
plt.imshow(galaxy)
plt.show()

#############################################################################
# Generate noisy observation
# --------------------------
# %%
# Generate a noisy observation
# ----------------------------
#
# Convolve the image with a point spread function (PSF) using the `convolve`
# function. Then add random Gaussian noise with standard deviation 0.0005
# using the `add_noise` function.
# To simulate an observation we convolve the image with a point spread function
# (PSF) using the ModOpt :py:func:`convolve <modopt.math.convolve.convolve>`
# function. Then we add random Gaussian noise with standard deviation
# ``0.0005`` using the ModOpt
# :py:func:`add_noise <modopt.signal.noise.add_noise>` function.

obs_data = add_noise(convolve(galaxy.data, psf.data), sigma=0.0005)

#############################################################################
# Create a PySAP image object
# %%
# Now we can show the noisy and blurred galaxy image.

image_obs = Image(data=np.abs(obs_data))
plt.imshow(obs_data)
plt.show()

#############################################################################
# Show the noisy galaxy image

plt.imshow(image_obs)

#############################################################################
# %%
# Deconvolve
# ----------
#
# Use the `sparse_deconv_condatvu` function to deconvolve the noisy image and
# set the maximum number of iterations to 3000.
# We will use the :py:func:`sparse_deconv_condatvu
# <astro.deconvolution.deconvolve.sparse_deconv_condatvu>` function from
# PySAP-Astro to deconvolve the noisy image. We set the maximum number of
# iterations for this function to ``3000``.

deconv_data = sparse_deconv_condatvu(obs_data, psf.data, n_iter=3000)

#############################################################################
# Create a PySAP image object for the result

image_rec = Image(data=np.abs(deconv_data))
# %%
# We can show the deconvolved galaxy image.

#############################################################################
# Show the deconvolved galaxy image
plt.imshow(deconv_data)
plt.show()

plt.imshow(image_rec)

#############################################################################
# %%
# Residual
# --------
#
# Create a PySAP image object for the residual
# Next, we calculate the residual of our deconvolved image to get a measure of
# how well the deconvolution process performed.

residual = Image(data=np.abs(galaxy.data - deconv_data))
residual = np.abs(galaxy.data - deconv_data)

#############################################################################
# Show the residual
# %%
# Finally, we can show the residual.

plt.imshow(residual)
plt.show()

# %%
# .. tip::
#
# Typically for a denoising problem we are aiming for a residual without any
# structure, i.e. just noise.
88 changes: 49 additions & 39 deletions examples/plot_galaxy_denoising.py
Original file line number Diff line number Diff line change
@@ -1,85 +1,95 @@
# -*- coding: utf-8 -*-
"""
Galaxy Image Denoising
======================
In this tutorial we will remove the noise from an example galaxy image.
.. codeauthor:: Samuel Farrens <[email protected]>
Import Dependencies
In this tutorial we will remove the noise from an example galaxy image using
the PySAP-Astro plug-in.
Import dependencies
-------------------
Import functions from PySAP and ModOpt.
Import functions from PySAP and
`ModOpt <https://cea-cosmic.github.io/ModOpt/>`_.
"""

import numpy as np
import matplotlib.pyplot as plt
from pysap import Image
from pysap.data import get_sample_data
from astro.denoising.denoise import denoise
from modopt.signal.noise import add_noise

#############################################################################
# Load the image of galaxy NGC2997
# %%
# Clean image
# -----------
#
# First, we load the image of galaxy NGC2997 from the PySAP sample data sets.

galaxy = get_sample_data('astro-ngc2997')

#############################################################################
# Show the clean galaxy image
# %%
# Then we can show the clean galaxy image that we will attempt to recover.

plt.imshow(galaxy)
plt.show()

#############################################################################
# Generate noisy observation
# --------------------------
# %%
# Generate a noisy observation
# ----------------------------
#
# Add random Gaussian noise with standard deviation 100 using the `add_noise`
# function.
# To simulate an observation of NGC2997 we add random Gaussian noise with
# standard deviation ``100`` using the
# ModOpt :py:func:`add_noise <modopt.signal.noise.add_noise>` function.

obs_data = add_noise(galaxy.data, sigma=100)

#############################################################################
# Create a PySAP image object

image_obs = Image(data=np.abs(obs_data))
# %%
# Now we can show the noisy galaxy image.

#############################################################################
# Show the noisy galaxy image

plt.imshow(image_obs)
plt.imshow(obs_data)
plt.show()

#############################################################################
# %%
# .. note::
#
# :math:`\sigma=100` is quite excesive, but we can more easily visualise the
# noise added to the image in this example.

# %%
# Denoise
# -------
#
# Use the `sparse_deconv_condatvu` function to deconvolve the noisy image and
# set the maximum number of iterations to 3000.

denoise_data = denoise(obs_data, n_scales=4)

#############################################################################
# Create a PySAP image object for the result
# We will use the :py:func:`denoise <astro.denoising.denoise.denoise>` function
# from PySAP-Astro to denoise the noisy image. We set the number of wavelet
# scales to ``4``.

image_rec = Image(data=np.abs(denoise_data))
denoised_data = denoise(obs_data, n_scales=4)

#############################################################################
# Show the deconvolved galaxy image
# %%
# We can show the denoised galaxy image.

plt.imshow(image_rec)
plt.imshow(denoised_data)
plt.show()

#############################################################################
# %%
# Residual
# --------
#
# Create a PySAP image object for the residual
# Next, we calculate the residual of our denoised image to get a measure of
# how well the denoising process performed.

residual = Image(data=np.abs(galaxy.data - denoise_data))
residual = np.abs(galaxy.data - denoised_data)

#############################################################################
# Show the residual
# %%
# Finally, we can show the residual.

plt.imshow(residual)
plt.show()

# %%
# .. tip::
#
# Typically for a denoising problem we are aiming for a residual without any
# structure, i.e. just noise.

0 comments on commit 5f75837

Please sign in to comment.