diff --git a/CHANGES.rst b/CHANGES.rst index 4ce7192f..430d7354 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -7,6 +7,13 @@ Bug Fixes - Improved errors/warnings when background region extends beyond bounds of image [#127] - Fixed boxcar weighting bug that often resulted in peak pixels having weight above 1 and erroneously triggered overlapping background errors [#125] +- Fixed boxcar weighting to handle zero width and edge of image cases [#141] + +New Features +^^^^^^^^^^^^ +- ``Background`` has new methods for exposing the 1D spectrum of the background or + background-subtracted regions [#143] + 1.1.0 ----- @@ -15,8 +22,6 @@ New Features ^^^^^^^^^^^^ - ``peak_method`` as an optional argument to ``KosmosTrace`` [#115] -- ``Background`` has new methods for exposing the 1D spectrum of the background or - background-subtracted regions [#143] API Changes ^^^^^^^^^^^ diff --git a/specreduce/background.py b/specreduce/background.py index 7345a039..47ccf9f4 100644 --- a/specreduce/background.py +++ b/specreduce/background.py @@ -86,15 +86,25 @@ def _to_trace(trace): raise ValueError('trace_object.trace_pos must be >= 1') return trace + if self.width < 0: + raise ValueError("width must be positive") + + if self.width == 0: + self.bkg_array = np.zeros(self.image.shape[self.disp_axis]) + return + bkg_wimage = np.zeros_like(self.image, dtype=np.float64) for trace in self.traces: trace = _to_trace(trace) - if (np.any(trace.trace.data >= self.image.shape[self.crossdisp_axis]) or - np.any(trace.trace.data < 0)): - raise ValueError("center of background window goes beyond image boundaries") - elif (np.any(trace.trace.data + self.width/2. >= self.image.shape[self.crossdisp_axis]) - or np.any(trace.trace.data - self.width/2. < 0)): - warnings.warn("background window extends beyond image boundaries") + 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]: + warnings.warn("background window extends beyond image boundaries " + + f"({windows_max} >= {self.image.shape[self.crossdisp_axis]})") + if windows_min < 0: + warnings.warn("background window extends beyond image boundaries " + + f"({windows_min} < 0)") + # pass trace.trace.data to ignore any mask on the trace bkg_wimage += _ap_weight_image(trace, self.width, @@ -104,18 +114,22 @@ def _to_trace(trace): if np.any(bkg_wimage > 1): raise ValueError("background regions overlapped") + if np.any(np.sum(bkg_wimage, axis=self.crossdisp_axis) == 0): + raise ValueError("background window does not remain in bounds across entire dispersion axis") # noqa if self.statistic == 'median': # make it clear in the expose image that partial pixels are fully-weighted bkg_wimage[bkg_wimage > 0] = 1 self.bkg_wimage = bkg_wimage + if self.statistic == 'average': - self.bkg_array = np.average(self.image, weights=self.bkg_wimage, axis=0) + self.bkg_array = np.average(self.image, weights=self.bkg_wimage, + axis=self.crossdisp_axis) elif self.statistic == 'median': med_image = self.image.copy() med_image[np.where(self.bkg_wimage) == 0] = np.nan - self.bkg_array = np.nanmedian(med_image, axis=0) + self.bkg_array = np.nanmedian(med_image, axis=self.crossdisp_axis) else: raise ValueError("statistic must be 'average' or 'median'") diff --git a/specreduce/extract.py b/specreduce/extract.py index 56d2b211..9965794f 100644 --- a/specreduce/extract.py +++ b/specreduce/extract.py @@ -43,22 +43,37 @@ def _get_boxcar_weights(center, hwidth, npix): A 2D image with weights assigned to pixels that fall within the defined aperture. """ - weights = np.zeros(npix) - - # shift center from integer to pixel space, where pixel N is [N-0.5, N+0.5), - # not [N, N+1). a pixel's integer index corresponds to its middle, not edge - center += 0.5 - - # pixels given full weight because they sit entirely within the aperture - fullpixels = [max(0, int(np.ceil(center - hwidth))), - min(int(np.floor(center + hwidth)), npix)] - weights[fullpixels[0]:fullpixels[1]] = 1 - - # pixels at the edges of the boxcar with partial weight, if any - if fullpixels[0] > 0: - weights[fullpixels[0] - 1] = hwidth - (center - fullpixels[0]) - if fullpixels[1] < npix: - weights[fullpixels[1]] = hwidth - (fullpixels[1] - center) + weights = np.zeros((npix)) + if hwidth == 0: + # the logic below would return all zeros anyways, so might as well save the time + # (negative widths should be avoided by earlier logic!) + return weights + + if center-hwidth > npix-0.5 or center+hwidth < -0.5: + # entire window is out-of-bounds + return weights + + lower_edge = max(-0.5, center-hwidth) # where -0.5 is lower bound of the image + upper_edge = min(center+hwidth, npix-0.5) # where npix-0.5 is upper bound of the image + + # let's avoid recomputing the round repeatedly + int_round_lower_edge = int(round(lower_edge)) + int_round_upper_edge = int(round(upper_edge)) + + # inner pixels that get full weight + # the round in conjunction with the +1 handles the half-pixel "offset", + # the upper bound doesn't have the +1 because array slicing is inclusive on the lower index and + # exclusive on the upper-index + # NOTE: round(-0.5) == 0, which is helpful here for the case where lower_edge == -0.5 + weights[int_round_lower_edge+1:int_round_upper_edge] = 1 + + # handle edge pixels (for cases where an edge pixel is fully-weighted, this will set it again, + # but should still compute a weight of 1. By using N:N+1, we avoid index errors if the edge + # is outside the image bounds. But we do need to avoid negative indices which would count + # from the end of the array. + if int_round_lower_edge >= 0: + weights[int_round_lower_edge:int_round_lower_edge+1] = round(lower_edge) + 0.5 - lower_edge + weights[int_round_upper_edge:int_round_upper_edge+1] = upper_edge - (round(upper_edge) - 0.5) return weights @@ -185,6 +200,9 @@ def __call__(self, image=None, trace_object=None, width=None, if trace_object.trace_pos < 1: raise ValueError('trace_object.trace_pos must be >= 1') + if width < 0: + raise ValueError("width must be positive") + # weight image to use for extraction wimage = _ap_weight_image( trace_object, diff --git a/specreduce/tests/test_background.py b/specreduce/tests/test_background.py index f5808fb1..28e4c3a7 100644 --- a/specreduce/tests/test_background.py +++ b/specreduce/tests/test_background.py @@ -6,7 +6,7 @@ from specutils import Spectrum1D from specreduce.background import Background -from specreduce.tracing import FlatTrace +from specreduce.tracing import FlatTrace, ArrayTrace # NOTE: same test image as in test_extract.py @@ -49,13 +49,36 @@ def test_background(): assert isinstance(bkg_spec, Spectrum1D) sub_spec = bg1.sub_spectrum() assert isinstance(sub_spec, Spectrum1D) + # test that width==0 results in no background + bg = Background.two_sided(image, trace, bkg_sep, width=0) + assert np.all(bg.bkg_image() == 0) -def test_oob(): +def test_warnings_errors(): # image.shape (30, 10) with pytest.warns(match="background window extends beyond image boundaries"): Background.two_sided(image, 25, 4, width=3) - with pytest.raises(ValueError, - match="center of background window goes beyond image boundaries"): - Background.two_sided(image, 25, 6, width=3) + # bottom of top window near/on top-edge of image (these should warn, but not fail) + with pytest.warns(match="background window extends beyond image boundaries"): + Background.two_sided(image, 25, 8, width=5) + + with pytest.warns(match="background window extends beyond image boundaries"): + Background.two_sided(image, 25, 8, width=6) + + with pytest.warns(match="background window extends beyond image boundaries"): + Background.two_sided(image, 25, 8, width=7) + + 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 + 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 + # 20 + 10 - 3 = 27 (lower edge of window on-image at right side of trace) + # 29 + 10 - 3 = 36 (lower edge of window off-image at right side of trace) + Background.one_sided(image, trace, 10, width=3) + + with pytest.raises(ValueError, match="width must be positive"): + Background.two_sided(image, 25, 2, width=-1)