From 787dd870dbcb3740edf71448a3c94ff944d89af3 Mon Sep 17 00:00:00 2001 From: AJANI Virginia Date: Fri, 11 Jun 2021 15:30:44 +0200 Subject: [PATCH 1/3] add script starlet l1-norm --- lenspack/starlet_l1norm.py | 89 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 89 insertions(+) create mode 100644 lenspack/starlet_l1norm.py diff --git a/lenspack/starlet_l1norm.py b/lenspack/starlet_l1norm.py new file mode 100644 index 0000000..3e5da6d --- /dev/null +++ b/lenspack/starlet_l1norm.py @@ -0,0 +1,89 @@ +# -*- coding: utf-8 -*- + +"""STARLET L1-NORM MODULE + +This module contains functions for computing the starlet l1norm +as defined in Eq. (1) of https://arxiv.org/pdf/2101.01542.pdf. + +""" + +import numpy as np +from astropy.stats import mad_std +from lenspack.image.transforms import starlet2d + + +def noise_coeff(image, nscales): + + """Compute the noise coefficients $\sigma_j$ + to get the estimate of the noise at the scale j + following Starck and Murtagh (1998). + + Parameters + ---------- + image : array_like + Two-dimensional input image. + nscales : int + Number of wavelet scales to compute. Should not exceed log2(N), where + N is the smaller of the two input dimensions. + + Returns + ------- + sigma_j : numpy.ndarray + Values of the standard deviation of the noise at scale j + """ + + noise_sigma=np.random.randn(image.shape[0],image.shape[0]) + noise_wavelet=starlet2d(noise_sigma, nscales) + sigma_j=np.array([np.std(scale) for scale in noise_wavelet]) + return sigma_j + + +def get_l1norm_noisy(image, noise, nscales, nbins): + + """Compute the starlet $\ell_1$-norm of a noisy image + following Eq. (1) of https://arxiv.org/abs/2101.01542. + + Parameters + ---------- + image : array_like + Two-dimensional input noiseless image. + noise : array_like + Two-dimensional input of the noise to be added to image + nscales : int + Number of wavelet scales to compute. Should not exceed log2(N), where + N is the smaller of the two input dimensions. + nbins : int + Number of bins in S/N desired for the summary statistic + + Returns + ------- + bins_snr, starlet_l1norm : tuple of 1D numpy arrays + Bin centers in S/N and Starlet $\ell_1$-norm of the noisy image + """ + + #add noise to noiseless image + image_noisy=image+noise + + #perform starlet decomposition + image_starlet = starlet2d(image_noisy, nscales) + + # estimate of the noise + noise_estimate = mad_std(image_noisy) + + std_coeff=noise_coeff(image, nscales) + + l1_coll = [] + bins_coll = [] + for image_temp, std_co in zip(image_starlet, std_coeff): + + std_scalej = std_co*noise_estimate + + snr = image_temp/std_scalej + thresholds_snr = np.linspace(np.min(snr), np.max(snr), nbins+1) + bins_snr=0.5*(thresholds_snr[:-1] + thresholds_snr[1:]) + digitized = np.digitize(snr, thresholds_snr) + bin_l1_norm = [np.sum(np.abs(snr[digitized == i])) for i in range(1, len(thresholds_snr))] + l1_coll.append(bin_l1_norm) + bins_coll.append(bins_snr) + + return np.array(bins_coll), np.array(l1_coll) From c2f6cedead8d93419ecf40c26355fae8aea3c355 Mon Sep 17 00:00:00 2001 From: AJANI Virginia Date: Fri, 11 Jun 2021 17:08:16 +0200 Subject: [PATCH 2/3] solve PEP8 errors --- lenspack/starlet_l1norm.py | 40 ++++++++++++++------------------------ 1 file changed, 15 insertions(+), 25 deletions(-) diff --git a/lenspack/starlet_l1norm.py b/lenspack/starlet_l1norm.py index 3e5da6d..2bd31e1 100644 --- a/lenspack/starlet_l1norm.py +++ b/lenspack/starlet_l1norm.py @@ -13,11 +13,9 @@ def noise_coeff(image, nscales): - - """Compute the noise coefficients $\sigma_j$ + """Compute the noise coefficients :math:`\sigma_j` to get the estimate of the noise at the scale j following Starck and Murtagh (1998). - Parameters ---------- image : array_like @@ -25,24 +23,20 @@ def noise_coeff(image, nscales): nscales : int Number of wavelet scales to compute. Should not exceed log2(N), where N is the smaller of the two input dimensions. - Returns ------- sigma_j : numpy.ndarray Values of the standard deviation of the noise at scale j """ - - noise_sigma=np.random.randn(image.shape[0],image.shape[0]) - noise_wavelet=starlet2d(noise_sigma, nscales) - sigma_j=np.array([np.std(scale) for scale in noise_wavelet]) + noise_sigma = np.random.randn(image.shape[0], image.shape[0]) + noise_wavelet = starlet2d(noise_sigma, nscales) + sigma_j = np.array([np.std(scale) for scale in noise_wavelet]) return sigma_j def get_l1norm_noisy(image, noise, nscales, nbins): - """Compute the starlet $\ell_1$-norm of a noisy image - following Eq. (1) of https://arxiv.org/abs/2101.01542. - + following Eq. (1) of https://arxiv.org/abs/2101.01542. Parameters ---------- image : array_like @@ -54,36 +48,32 @@ def get_l1norm_noisy(image, noise, nscales, nbins): N is the smaller of the two input dimensions. nbins : int Number of bins in S/N desired for the summary statistic - Returns ------- bins_snr, starlet_l1norm : tuple of 1D numpy arrays Bin centers in S/N and Starlet $\ell_1$-norm of the noisy image """ - #add noise to noiseless image - image_noisy=image+noise - - #perform starlet decomposition + # add noise to noiseless image + image_noisy = image + noise + # perform starlet decomposition image_starlet = starlet2d(image_noisy, nscales) - # estimate of the noise noise_estimate = mad_std(image_noisy) - - std_coeff=noise_coeff(image, nscales) + std_coeff = noise_coeff(image, nscales) l1_coll = [] bins_coll = [] for image_temp, std_co in zip(image_starlet, std_coeff): - std_scalej = std_co*noise_estimate + std_scalej = std_co * noise_estimate - snr = image_temp/std_scalej - thresholds_snr = np.linspace(np.min(snr), np.max(snr), nbins+1) - bins_snr=0.5*(thresholds_snr[:-1] + thresholds_snr[1:]) + snr = image_temp / std_scalej + thresholds_snr = np.linspace(np.min(snr), np.max(snr), nbins + 1) + bins_snr = 0.5 * (thresholds_snr[:-1] + thresholds_snr[1:]) digitized = np.digitize(snr, thresholds_snr) - bin_l1_norm = [np.sum(np.abs(snr[digitized == i])) for i in range(1, len(thresholds_snr))] + bin_l1_norm = [np.sum(np.abs(snr[digitized == i])) + for i in range(1, len(thresholds_snr))] l1_coll.append(bin_l1_norm) bins_coll.append(bins_snr) - return np.array(bins_coll), np.array(l1_coll) From 3f18677acce73e168032541c61d18d28539a93cb Mon Sep 17 00:00:00 2001 From: AJANI Virginia Date: Mon, 14 Jun 2021 15:37:44 +0200 Subject: [PATCH 3/3] change names variables after first review --- lenspack/starlet_l1norm.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/lenspack/starlet_l1norm.py b/lenspack/starlet_l1norm.py index 2bd31e1..9efcb23 100644 --- a/lenspack/starlet_l1norm.py +++ b/lenspack/starlet_l1norm.py @@ -13,7 +13,7 @@ def noise_coeff(image, nscales): - """Compute the noise coefficients :math:`\sigma_j` + """Compute the noise coefficients :math:`\sigma^{e}_{j}` to get the estimate of the noise at the scale j following Starck and Murtagh (1998). Parameters @@ -25,17 +25,17 @@ def noise_coeff(image, nscales): N is the smaller of the two input dimensions. Returns ------- - sigma_j : numpy.ndarray + coeff_j : numpy.ndarray Values of the standard deviation of the noise at scale j """ noise_sigma = np.random.randn(image.shape[0], image.shape[0]) noise_wavelet = starlet2d(noise_sigma, nscales) - sigma_j = np.array([np.std(scale) for scale in noise_wavelet]) - return sigma_j + coeff_j = np.array([np.std(scale) for scale in noise_wavelet]) + return coeff_j def get_l1norm_noisy(image, noise, nscales, nbins): - """Compute the starlet $\ell_1$-norm of a noisy image + """Compute the starlet :math:`\ell_1`-norm of a noisy image following Eq. (1) of https://arxiv.org/abs/2101.01542. Parameters ---------- @@ -51,7 +51,7 @@ def get_l1norm_noisy(image, noise, nscales, nbins): Returns ------- bins_snr, starlet_l1norm : tuple of 1D numpy arrays - Bin centers in S/N and Starlet $\ell_1$-norm of the noisy image + Bin centers in S/N and Starlet :math:`\ell_1`-norm of the noisy image """ # add noise to noiseless image @@ -60,15 +60,15 @@ def get_l1norm_noisy(image, noise, nscales, nbins): image_starlet = starlet2d(image_noisy, nscales) # estimate of the noise noise_estimate = mad_std(image_noisy) - std_coeff = noise_coeff(image, nscales) + coeff_j = noise_coeff(image, nscales) l1_coll = [] bins_coll = [] - for image_temp, std_co in zip(image_starlet, std_coeff): + for image_j, std_co in zip(image_starlet, coeff_j): - std_scalej = std_co * noise_estimate + sigma_j = std_co * noise_estimate - snr = image_temp / std_scalej + snr = image_j / sigma_j thresholds_snr = np.linspace(np.min(snr), np.max(snr), nbins + 1) bins_snr = 0.5 * (thresholds_snr[:-1] + thresholds_snr[1:]) digitized = np.digitize(snr, thresholds_snr)