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.

Sign up for GitHub

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

Jump to bottom

Update utils.py - add Median Absolute Deviation option to radial_profile #591

Merged
merged 2 commits into from
Dec 11, 2023
Merged

Update utils.py - add Median Absolute Deviation option to radial_profile #591

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
78b2537
Update utils.py
saikrishanthpm Oct 5, 2023
7b1babb
maintainer edits for PR #591: add unit test, some docs clarifications
mperrin Dec 11, 2023
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
9 changes: 9 additions & 0 deletions poppy/tests/test_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,15 @@ def test_radial_profile(plot=False):
assert np.allclose(prof2, prof3)
# TODO compare those to be near a sinc profile as expected?

# Test computing the deviation profiles
rad4, prof_stds = poppy.radial_profile(psf, stddev=True)
rad4, prof_mads = poppy.radial_profile(psf, mad=True)
# TODO add some more rigorous test of correctness; for now just test
# dimensionality of the outputs
assert prof_stds.shape == prof.shape, "Radial profile and stddev array output sizes should match"
assert prof_mads.shape == prof.shape, "Radial profile and median absolute deviation array output sizes should match"


def test_radial_profile_of_offset_source():
"""Test that we can compute the radial profile for a source slightly outside the FOV,
compare that to a calculation for a centered source, and check we get consistent results
Expand Down
28 changes: 24 additions & 4 deletions poppy/utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -539,8 +539,8 @@ def display_profiles(hdulist_or_filename=None, ext=0, overplot=False, title=None
plt.text(ee_lev + 0.1, level + yoffset, 'EE=%2d%% at r=%.3f"' % (level * 100, ee_lev))


def radial_profile(hdulist_or_filename=None, ext=0, ee=False, center=None, stddev=False, binsize=None, maxradius=None,
normalize='None', pa_range=None, slice=0):
def radial_profile(hdulist_or_filename=None, ext=0, ee=False, center=None, stddev=False, mad=False,
binsize=None, maxradius=None, normalize='None', pa_range=None, slice=0):
""" Compute a radial profile of the image.

This computes a discrete radial profile evaluated on the provided binsize. For a version
Expand All @@ -563,6 +563,9 @@ def radial_profile(hdulist_or_filename=None, ext=0, ee=False, center=None, stdde
size of step for profile. Default is pixel size.
stddev : bool
Compute standard deviation in each radial bin, not average?
mad : bool
Compute median absolute deviation (MAD) in each radial bin.
Cannot be used at same time as stddev; pick one or the other.
normalize : string
set to 'peak' to normalize peak intensity =1, or to 'total' to normalize total flux=1.
Default is no normalization (i.e. retain whatever normalization was used in computing the PSF itself)
Expand All @@ -583,7 +586,9 @@ def radial_profile(hdulist_or_filename=None, ext=0, ee=False, center=None, stdde
Tuple containing (radius, profile) or (radius, profile, EE) depending on what is requested.
The radius gives the center radius of each bin, while the EE is given inside the whole bin
so you should use (radius+binsize/2) for the radius of the EE curve if you want to be
as precise as possible.
as precise as possible. The profile will be either the average within each bin, or the
standard or median absolute deviation within each bin if one of those options is selected.

"""
if isinstance(hdulist_or_filename, str):
hdu_list = fits.open(hdulist_or_filename)
Expand Down Expand Up @@ -668,6 +673,7 @@ def radial_profile(hdulist_or_filename=None, ext=0, ee=False, center=None, stdde
radialprofile2 = radialprofile2[crop]

if stddev:
# Compute standard deviation in each radial bin
stddevs = np.zeros_like(radialprofile2)
r_pix = r * binsize
for i, radius in enumerate(rr):
Expand All @@ -679,9 +685,23 @@ def radial_profile(hdulist_or_filename=None, ext=0, ee=False, center=None, stdde
stddevs[i] = np.nanstd(image[wg])
return rr, stddevs

if not ee:
elif mad:
# Compute median absolute deviation in each radial bin
mads = np.zeros_like(radialprofile2)
r_pix = r * binsize
for i, radius in enumerate(rr):
if i == 0:
wg = np.where(r < radius + binsize / 2)
else:
wg = np.where((r_pix >= (radius - binsize / 2)) & (r_pix < (radius + binsize / 2)))
mads[i] = np.nanmedian(np.absolute(image[wg]-np.nanmedian(image[wg])))
return rr, mads

elif not ee:
# (Default behavior) Compute average in each radial bin
return rr, radialprofile2
else:
# also return the cumulative sum within each radial bin, i.e. the encircled energy
ee = csim[rind]
return rr, radialprofile2, ee

Expand Down