Skip to content

Commit

Permalink
.
Browse files Browse the repository at this point in the history
  • Loading branch information
cshanahan1 committed Apr 22, 2024
1 parent c16abd4 commit 494f618
Show file tree
Hide file tree
Showing 6 changed files with 384 additions and 188 deletions.
25 changes: 20 additions & 5 deletions specreduce/background.py
Original file line number Diff line number Diff line change
Expand Up @@ -98,16 +98,27 @@ def __post_init__(self):
dispersion axis
crossdisp_axis : int
cross-dispersion axis
mask_treatment : string
The method for handling masked or non-finite data. Choice of `filter`,
`omit`, or `zero-fill`. If `filter` is chosen, masked/non-finite data
will be filtered during the fit to each bin/column (along disp. axis) to
find the peak. If `omit` is chosen, columns along disp_axis with any
masked/non-finite data values will be fully masked (i.e, 2D mask is
collapsed to 1D and applied). If `zero-fill` is chosen, masked/non-finite
data will be replaced with 0.0 in the input image, and the mask will then
be dropped. For all three options, the input mask (optional on input
NDData object) will be combined with a mask generated from any non-finite
values in the image data.
"""

# Parse image, including masked/nonfinite data handling based on
# choice of `mask_treatment`. returns a Spectrum1D
# choice of `mask_treatment`. Any uncaught nonfinte data values will be
# masked as well. Returns a Spectrum1D.
self.image = self._parse_image(self.image)

# _parse_image returns a Spectrum1D. convert this to a masked array
# for ease of calculations here (even if there is no masked data).
# Note: uncertainties are dropped, this should also be addressed at
# some point probably across the package.
# always work with masked array, even if there is no masked
# or nonfinite data, in case padding is needed. if not, mask will be
# dropped at the end and a regular array will be returned.
img = np.ma.masked_array(self.image.data, self.image.mask)

if self.width < 0:
Expand All @@ -120,6 +131,9 @@ def __post_init__(self):

bkg_wimage = np.zeros_like(self.image.data, dtype=np.float64)
for trace in self.traces:
# note: ArrayTrace can have masked values, but if it does a MaskedArray
# will be returned so this should be reflected in the window size here
# (i.e, np.nanmax is not required.)
windows_max = trace.trace.data.max() + self.width/2
windows_min = trace.trace.data.min() - self.width/2
if windows_max >= self.image.shape[self.crossdisp_axis]:
Expand All @@ -136,6 +150,7 @@ def __post_init__(self):
self.crossdisp_axis,
self.image.shape)


if np.any(bkg_wimage > 1):
raise ValueError("background regions overlapped")
if np.any(np.sum(bkg_wimage, axis=self.crossdisp_axis) == 0):
Expand Down
1 change: 0 additions & 1 deletion specreduce/core.py
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,6 @@ def _mask_and_nonfinite_data_handling(self, image, mask):
# nothing needs to be done. input mask (combined with nonfinite data)
# remains with data as-is.


if mask_treatment == 'zero-fill':
# make a copy of the input image since we will be modifying it
image = deepcopy(image)
Expand Down
83 changes: 44 additions & 39 deletions specreduce/tests/test_background.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,13 +24,13 @@ def test_background(mk_test_img_raw, mk_test_spec_no_spectral_axis,

# all the following should be equivalent, whether image's spectral axis
# is in pixels or physical units:
bg1 = Background(image, [trace-bkg_sep, trace+bkg_sep], width=bkg_width)
bg1 = Background(image, [trace - bkg_sep, trace + bkg_sep], width=bkg_width)
bg2 = Background.two_sided(image, trace, bkg_sep, width=bkg_width)
bg3 = Background.two_sided(image, trace_pos, bkg_sep, width=bkg_width)
assert np.allclose(bg1.bkg_image().flux, bg2.bkg_image().flux)
assert np.allclose(bg1.bkg_image().flux, bg3.bkg_image().flux)

bg4 = Background(image_um, [trace-bkg_sep, trace+bkg_sep], width=bkg_width)
bg4 = Background(image_um, [trace - bkg_sep, trace + bkg_sep], width=bkg_width)
bg5 = Background.two_sided(image_um, trace, bkg_sep, width=bkg_width)
bg6 = Background.two_sided(image_um, trace_pos, bkg_sep, width=bkg_width)
assert np.allclose(bg1.bkg_image().flux, bg4.bkg_image().flux)
Expand Down Expand Up @@ -74,7 +74,7 @@ def test_background(mk_test_img_raw, mk_test_spec_no_spectral_axis,
stats = ['average', 'median']

for st in stats:
bg = Background(img, trace-bkg_sep, width=bkg_width, statistic=st)
bg = Background(img, trace - bkg_sep, width=bkg_width, statistic=st)
assert np.isnan(bg.image.flux).sum() == 2
assert np.isnan(bg._bkg_array).sum() == 0
assert np.isnan(bg.bkg_spectrum().flux).sum() == 0
Expand All @@ -101,7 +101,7 @@ def test_warnings_errors(mk_test_spec_no_spectral_axis):
with pytest.warns(match="background window extends beyond image boundaries"):
Background.two_sided(image, 7, 5, width=6)

trace = ArrayTrace(image, trace=np.arange(10)+20) # from 20 to 29
trace = ArrayTrace(image, trace=np.arange(10) + 20) # from 20 to 29
with pytest.warns(match="background window extends beyond image boundaries"):
with pytest.raises(ValueError,
match="background window does not remain in bounds across entire dispersion axis"): # noqa
Expand All @@ -112,8 +112,8 @@ def test_warnings_errors(mk_test_spec_no_spectral_axis):
with pytest.raises(ValueError, match="width must be positive"):
Background.two_sided(image, 25, 2, width=-1)

def test_trace_inputs(mk_test_img_raw):

def test_trace_inputs(mk_test_img_raw):
"""
Tests for the input argument 'traces' to `Background`. This should accept
a list of or a single Trace object, or a list of or a single (positive)
Expand Down Expand Up @@ -143,78 +143,84 @@ def test_trace_inputs(mk_test_img_raw):
with pytest.raises(ValueError, match=match_str):
Background(image, 'non_valid_trace_pos')


class TestMasksBackground():

"""
Various test functions to test how masked and non-finite data is handled
in `Background.
in `Background. There are three currently implemented options for masking
in Background: filter, omit, and zero-fill.
"""

def mk_img(self, nrows=4, ncols=5, nan_slices=None):
"""
Make a simpleimage to test masking in Background.
Optionally add NaNs to data. Returned array is in u.DN.
Make a simple gradient image to test masking in Background.
Optionally add NaNs to data with `nan_slices`. Returned array is in
u.DN.
"""

img = np.tile((np.arange(1., ncols+1)), (nrows, 1))
img = np.tile((np.arange(1., ncols + 1)), (nrows, 1))

if nan_slices: # add nans in data
for s in nan_slices:
img[s] = np.nan

return img * u.DN

def test_fully_masked_column(self):
@pytest.mark.parametrize("mask", ["filter", "omit", "zero-fill"])
def test_fully_masked_column(self, mask):
"""
Test what happens when a full column is masked, not the entire
image. In this case, the background value for that fully-masked
column should be 0.0, with no error or warning raised.
Test background with some fully-masked columns (not fully masked image).
In this case, the background value for that fully-masked column should
be 0.0, with no error or warning raised. This is the case for
mask_treatment=filter, omit, or zero-fill.
"""

img = np.ones((12, 12))
img = self.mk_img(nrows=10, ncols=10)
img[:, 0:1] = np.nan

bkg = Background(img, traces=FlatTrace(img, 6))

bkg = Background(img, traces=FlatTrace(img, 6), mask_treatment=mask)
assert np.all(bkg.bkg_image().data[:, 0:1] == 0.0)


def test_fully_masked(self):
@pytest.mark.parametrize("mask", ["filter", "omit", "zero-fill"])
def test_fully_masked_image(self, mask):
"""
Test that the appropriate error is raised by `Background` when image
is fully masked/NaN.
"""

with pytest.raises(ValueError, match='Image is fully masked.'):
# fully NaN image
img = np.zeros((4, 5)) * np.nan
Background(img, traces=FlatTrace(self.mk_img(), 2))
img = self.mk_img() * np.nan
Background(img, traces=FlatTrace(self.mk_img(), 2), mask_treatment=mask)

with pytest.raises(ValueError, match='Image is fully masked.'):
# fully masked image (should be equivilant)
# fully masked image (should be equivalent)
img = NDData(np.ones((4, 5)), mask=np.ones((4, 5)))
Background(img, traces=FlatTrace(self.mk_img(), 2))
Background(img, traces=FlatTrace(self.mk_img(), 2), mask_treatment=mask)

# Now test that an image that isn't fully masked, but is fully masked
# within the window determined by `width`, produces the correct result
# within the window determined by `width`, produces the correct result.
# only applicable for mask_treatment=filter, because this is the only
# option that allows a slice of masked values that don't span all rows.
msg = 'Image is fully masked within background window determined by `width`.'
with pytest.raises(ValueError, match=msg):
img = self.mk_img(nrows=12, ncols=12, nan_slices=[np.s_[3:10, :]])
Background(img, traces=FlatTrace(img, 6), width=7)

@pytest.mark.filterwarnings("ignore:background window extends beyond image boundaries")
@pytest.mark.parametrize("method,expected",
[("filter", np.array([1., 2., 3., 4., 5., 6., 7.,
[("filter", np.array([1., 2., 3., 4., 5., 6., 7.,
8., 9., 10., 11., 12.])),
("omit", np.array([0., 2., 3., 0., 5., 6.,
7., 0., 9., 10., 11., 12.])),
("zero-fill", np.array([ 0.58333333, 2., 3.,
2.33333333, 5., 6., 7.,
7.33333333, 9., 10., 11.,
12.]))])
("omit", np.array([0., 2., 3., 0., 5., 6.,
7., 0., 9., 10., 11., 12.])),
("zero-fill", np.array([0.58333333, 2., 3.,
2.33333333, 5., 6., 7.,
7.33333333, 9., 10., 11.,
12.]))])
def test_mask_treatment_bkg_img_spectrum(self, method, expected):
"""
This test function tests `Backgroud.bkg_image` and
"""
This test function tests `Backgroud.bkg_image` and
`Background.bkg_spectrum` when there is masked data. It also tests
background subtracting the image, and returning the spectrum of the
background subtracted image. This test is parameterized over all
Expand All @@ -227,8 +233,8 @@ def test_mask_treatment_bkg_img_spectrum(self, method, expected):

# make image, set some value to nan, which will be masked in the function
image1 = self.mk_img(nrows=img_size, ncols=img_size,
nan_slices=[np.s_[5:10, 0], np.s_[7:12, 3],
np.s_[2, 7]])
nan_slices=[np.s_[5:10, 0], np.s_[7:12, 3],
np.s_[2, 7]])

# also make an image that doesn't have nonf data values, but has
# masked values at the same locations, to make sure they give the same
Expand All @@ -240,13 +246,13 @@ def test_mask_treatment_bkg_img_spectrum(self, method, expected):
for image in [image1, image2]:

# construct a flat trace in center of image
trace = FlatTrace(image, img_size/2)
trace = FlatTrace(image, img_size / 2)

# create 'Background' object with `mask_treatment` set
# 'width' should be > size of image to use all pix (but warning will
# be raised, which we ignore.)
background = Background(image, mask_treatment=method,
traces=trace, width=img_size+1)
traces=trace, width=img_size + 1)

# test background image matches 'expected'
bk_img = background.bkg_image()
Expand All @@ -265,7 +271,7 @@ def test_sub_bkg_image(self):
"""
Test that masked and nonfinite data is handled correctly when subtracting
background from image, for all currently implemented masking
options ('filter', 'omit', and 'zero-fill').
options ('filter', 'omit', and 'zero-fill').
"""

# make image, set some value to nan, which will be masked in the function
Expand All @@ -288,7 +294,7 @@ def test_sub_bkg_image(self):
# 2d mask is reduced to a 1d mask to mask out full columns in the
# presence of any nans - this means that (as tested above in
# `test_mask_treatment_bkg_img_spectrum`) those columns will have 0.0
# background. In this case, image.mask is expanded to mask full
# background. In this case, image.mask is expanded to mask full
# columns - the image itself will not have full columns set to np.nan,
# so there are still valid background subtracted data values in this
# case, but the corresponding mask for that entire column will be masked.
Expand All @@ -314,4 +320,3 @@ def test_sub_bkg_image(self):

assert np.all(np.isfinite(subtracted_img_zero_fill.data))
assert np.all(subtracted_img_zero_fill.mask == 0)

60 changes: 60 additions & 0 deletions specreduce/tests/test_extract.py
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,7 @@ def test_boxcar_extraction(mk_test_img):
assert np.allclose(spectrum.flux.value, np.full_like(spectrum.flux.value, 67.15))



def test_boxcar_outside_image_condition(mk_test_img):
#
# Trace is such that extraction aperture lays partially outside the image
Expand Down Expand Up @@ -326,3 +327,62 @@ def test_horne_interpolated_nbins_fails(mk_test_img):
spatial_profile={'name': 'interpolated_profile',
'n_bins_interpolated_profile': 100})
ex.spectrum

class TestMasksExtract():

def mk_flat_gauss_img(nrows=200, ncols=160, nan_slices=None, add_noise=True):

"""
Makes a flat gaussian image for testing, with optional added gaussian
nosie and optional data values set to NaN. Variance is included, which
is required by HorneExtract. Returns a Spectrum1D with flux, spectral
axis, and uncertainty.
"""

sigma_pix = 4
col_model = models.Gaussian1D(amplitude=1, mean=nrows/2,
stddev=sigma_pix)
spec2dvar = np.ones((nrows, ncols))
noise = 0
if add_noise:
np.random.seed(7)
sigma_noise = 1
noise = np.random.normal(scale=sigma_noise, size=(nrows, ncols))

index_arr = np.tile(np.arange(nrows), (ncols, 1))
img = col_model(index_arr.T) + noise

if nan_slices: # add nans in data
for s in nan_slices:
img[s] = np.nan

wave = np.arange(0, img.shape[1], 1)
objectspec = Spectrum1D(spectral_axis=wave*u.m, flux=img*u.Jy,
uncertainty=VarianceUncertainty(spec2dvar*u.Jy*u.Jy))

return objectspec

def test_boxcar_fully_masked():
"""
Test that the appropriate error is raised by `BoxcarExtract` when image
is fully masked/NaN.
"""

img = mk_img()

with pytest.raises(ValueError, match='Image is fully masked.'):
# fully NaN image
img = np.zeros((4, 5)) * np.nan
Background(img, traces=FlatTrace(self.mk_img(), 2))

with pytest.raises(ValueError, match='Image is fully masked.'):
# fully masked image (should be equivalent)
img = NDData(np.ones((4, 5)), mask=np.ones((4, 5)))
Background(img, traces=FlatTrace(self.mk_img(), 2))

# Now test that an image that isn't fully masked, but is fully masked
# within the window determined by `width`, produces the correct result
msg = 'Image is fully masked within background window determined by `width`.'
with pytest.raises(ValueError, match=msg):
img = self.mk_img(nrows=12, ncols=12, nan_slices=[np.s_[3:10, :]])
Background(img, traces=FlatTrace(img, 6), width=7)
Loading

0 comments on commit 494f618

Please sign in to comment.