Rewrite boxcar weights logic (#141)

* rewrite boxcar weights to handle subpixel width and extending off image
* forbid negative widths
* loosen restriction on window extending off image: warning is still raised if any window extends off the image, error is only raised if the window is entirely off image at any point in the dispersion axis
* remove hardcoded axis=0 and reference crossdisp_axis instead
* handle case where bottom of window is within 1 pixel of top of image

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
 ^^^^^^^^^^^

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'")
 

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,

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)