TiltCorrection and TiltSolution improvements

CHANGES.rst

@@ -30,6 +30,13 @@ New Features
   (``AWAV-GRA``) and vacuum (``WAVE-GRI``) spectral axis types in FITS WCS
   export. [#316]
 
+- Added ``subtract_baseline`` and ``baseline_window`` options to
+  ``line_matching.find_arc_lines`` and ``TiltCorrection.find_arc_lines`` that estimate
+  the baseline flux of the spectrum with a sigma-clipped median (globally, or in running
+  windows) and remove it before the line detection. Baseline subtraction is on by default
+  in ``TiltCorrection.find_arc_lines``, which previously required background-subtracted
+  arc frames to detect any lines. [#XXX]
+
 API Changes
 ^^^^^^^^^^^
 
@@ -38,6 +45,21 @@ API Changes
   on instantiation, and its removal has been rescheduled from v2.0 to v1.11.
   Use ``specreduce.wavecal1d.WavelengthCalibration1D`` instead. [#316]
 
+Bug Fixes
+^^^^^^^^^
+
+- ``TiltCorrection.find_arc_lines`` now raises a ``ValueError`` when the reference row or
+  any cross-dispersion sample lies outside the image frame, and discards detected lines
+  whose fitted centroid is not finite or falls outside the frame along the dispersion
+  axis. Previously a sample past the frame raised an ``IndexError``, a negative sample
+  silently wrapped around to the opposite edge of the frame, and a failed line fit could
+  crash the nearest-neighbour matching. [#XXX]
+
+- ``TiltSolution.resample`` now propagates the input uncertainty, returned in the same
+  uncertainty class as the input, marks output bins that received a contribution from a
+  masked input pixel, and copies the input metadata to the resampled ``NDData``.
+  Previously all three were dropped. [#XXX]
+
 Other changes
 ^^^^^^^^^^^^^
 

docs/tilt_correction/tilt_correction.rst

@@ -118,8 +118,17 @@ Parameters:
 
 *   **fwhm**: Expected full width at half maximum of the spectral lines (in pixels), used by
     the line-finding algorithm.
-*   **noise_factor**: Multiplier for noise thresholding — lines below
-    ``noise_factor × noise_level`` are rejected. Default is 5.
+*   **noise_factor**: Multiplier for noise thresholding — pixels less than
+    ``noise_factor × uncertainty`` above the baseline are not considered line candidates.
+    Default is 5.
+*   **subtract_baseline**: Estimate the baseline (background) flux of each row with a
+    sigma-clipped median and remove it before the line detection. The detection compares
+    the flux against zero, so arc frames with a pedestal above ``noise_factor × uncertainty``
+    would otherwise yield no usable lines. Default is ``True``; set it to ``False`` for
+    frames that are already background-subtracted.
+*   **baseline_window**: Width in pixels of the chunks used for a running baseline
+    estimate, which lets the baseline follow a slowly varying background such as
+    scattered light. Default is ``None``, meaning a single global median per row.
 
 3. Fitting the Tilt Model
 *************************
@@ -262,4 +271,7 @@ flux-conserving resampling independently of the calibration workflow.
         corrected = ts.resample(science_frame, bin_edges=np.linspace(50, 950, 501))
 
     The ``resample`` method accepts a ``mask_treatment`` parameter with the same options as
-    the :class:`~specreduce.tilt_correction.TiltCorrection` constructor.
+    the :class:`~specreduce.tilt_correction.TiltCorrection` constructor. The returned
+    :class:`~astropy.nddata.NDData` carries the resampled uncertainty (in the same
+    uncertainty class as the input), a mask flagging every bin that received a contribution
+    from a masked input pixel, and a copy of the input metadata.

specreduce/line_matching.py

@@ -5,9 +5,10 @@
 import astropy.units as u
 import numpy as np
 from astropy.nddata import StdDevUncertainty
-from astropy.stats import gaussian_fwhm_to_sigma, gaussian_sigma_to_fwhm
+from astropy.stats import gaussian_fwhm_to_sigma, gaussian_sigma_to_fwhm, sigma_clipped_stats
 from astropy.modeling import models
 from astropy.table import QTable
+from astropy.utils.exceptions import AstropyUserWarning
 from astropy.wcs import WCS as astropy_WCS
 from gwcs.wcs import WCS as gWCS
 from specutils.fitting import find_lines_threshold, fit_lines
@@ -17,11 +18,85 @@
 __all__ = ["find_arc_lines", "match_lines_wcs"]
 
 
+def _estimate_baseline(
+    flux: np.ndarray, window: int | None = None, sigma: float = 3.0, maxiters: int = 5
+) -> np.ndarray:
+    """Estimate the baseline (pedestal) flux of an arc spectrum with a sigma-clipped median.
+
+    Emission lines are positive outliers, so the sigma-clipped median tracks the
+    line-free background level. When ``window`` is given, the median is computed
+    in consecutive chunks of ``window`` pixels and linearly interpolated between
+    the chunk centres, which lets the estimate follow a slowly varying background.
+
+    Parameters
+    ----------
+    flux
+        The spectrum flux values.
+    window
+        Approximate width of the chunks in pixels. If ``None``, or if the spectrum
+        is too short to hold at least two chunks, a single global median is used.
+    sigma
+        Clipping threshold in standard deviations passed to
+        `~astropy.stats.sigma_clipped_stats`.
+    maxiters
+        Maximum number of clipping iterations.
+
+    Returns
+    -------
+    numpy.ndarray
+        The baseline estimate, with the same shape as ``flux``.
+    """
+    flux = np.asarray(flux, dtype=float)
+    n = flux.size
+    global_median = sigma_clipped_stats(flux, sigma=sigma, maxiters=maxiters)[1]
+
+    if window is None:
+        return np.full(n, global_median)
+
+    if isinstance(window, bool) or int(window) != window or window <= 0:
+        raise ValueError("baseline window must be a positive integer number of pixels.")
+    window = int(window)
+
+    # Split the spectrum into equal-sized chunks of roughly ``window`` pixels.
+    nchunks = int(round(n / window))
+    if nchunks < 2:
+        return np.full(n, global_median)
+
+    edges = np.linspace(0, n, nchunks + 1).astype(int)
+    centres = 0.5 * (edges[:-1] + edges[1:] - 1)
+    medians = np.array(
+        [
+            sigma_clipped_stats(flux[lo:hi], sigma=sigma, maxiters=maxiters)[1]
+            for lo, hi in zip(edges[:-1], edges[1:])
+        ]
+    )
+    ok = np.isfinite(medians)
+    if ok.sum() < 2:
+        return np.full(n, global_median)
+    centres, medians = centres[ok], medians[ok]
+
+    # Interpolate between the chunk centres and extrapolate linearly beyond the
+    # outermost ones, so a smooth gradient is followed all the way to the edges.
+    x = np.arange(n)
+    baseline = np.interp(x, centres, medians)
+    left = x < centres[0]
+    right = x > centres[-1]
+    baseline[left] = medians[0] + (x[left] - centres[0]) * np.diff(medians[:2]) / np.diff(
+        centres[:2]
+    )
+    baseline[right] = medians[-1] + (x[right] - centres[-1]) * np.diff(medians[-2:]) / np.diff(
+        centres[-2:]
+    )
+    return baseline
+
+
 def find_arc_lines(
     spectrum: Spectrum,
     fwhm: float | u.Quantity = 5.0 * u.pix,
     window: float = 3.0,
     noise_factor: float = 5.0,
+    subtract_baseline: bool = True,
+    baseline_window: int | None = None,
 ) -> QTable:
     """
     Find arc lines in a spectrum using `~specutils.fitting.find_lines_threshold` and
@@ -30,7 +105,7 @@ def find_arc_lines(
     Parameters
     ----------
     spectrum : The extracted arc spectrum to search for lines. It should be background-subtracted
-        and must have an "uncertainty" attribute.
+        unless ``subtract_baseline`` is set, and must have an "uncertainty" attribute.
 
     fwhm
         Estimated full-width half-maximum of the lines in pixels.
@@ -42,10 +117,21 @@ def find_arc_lines(
         The factor to multiply the uncertainty by to determine the noise threshold
         in the `~specutils.fitting.find_lines_threshold` routine.
 
+    subtract_baseline
+        If ``True``, estimate the baseline (background) flux with a sigma-clipped median
+        and subtract it before the line detection. The threshold routine compares the flux
+        against zero, so a spectrum with a pedestal above ``noise_factor × uncertainty``
+        cannot be searched without this.
+
+    baseline_window
+        Width in pixels of the chunks used for a running baseline estimate. If ``None``,
+        a single global median is subtracted. Only used when ``subtract_baseline`` is ``True``.
+
     Returns
     -------
     QTable
         A table of detected arc lines and their properties: centroid, fwhm, and amplitude.
+        When the baseline is subtracted, the amplitudes are measured above the baseline.
     """
     # If fwhm is a float, convert it to a Quantity with the same unit as the spectral axis
     # of the input spectrum.
@@ -59,7 +145,23 @@ def find_arc_lines(
         spectrum = deepcopy(spectrum)
         spectrum.uncertainty = StdDevUncertainty(np.sqrt(np.abs(spectrum.flux.value)))
 
-    detected_lines = find_lines_threshold(spectrum, noise_factor=noise_factor)
+    if subtract_baseline:
+        baseline = _estimate_baseline(spectrum.flux.value, window=baseline_window)
+        spectrum = Spectrum(
+            (spectrum.flux.value - baseline) * spectrum.flux.unit,
+            spectral_axis=spectrum.spectral_axis,
+            uncertainty=spectrum.uncertainty,
+            mask=spectrum.mask,
+        )
+
+    with warnings.catch_warnings():
+        if subtract_baseline:
+            # The baseline is known to be removed, so specutils' continuum heuristic
+            # (median flux below 0.01 sigma) only adds noise for strong-line arc spectra.
+            warnings.filterwarnings(
+                "ignore", message="Spectrum is not below the threshold", category=AstropyUserWarning
+            )
+        detected_lines = find_lines_threshold(spectrum, noise_factor=noise_factor)
     detected_lines = detected_lines[detected_lines["line_type"] == "emission"]
 
     centroids = []

specreduce/tests/test_line_matching.py

@@ -5,12 +5,13 @@
 from astropy.wcs import WCS
 from astropy.modeling import models
 from astropy.nddata import StdDevUncertainty
+from astropy.utils.exceptions import AstropyUserWarning
 from specutils.fitting import fit_generic_continuum
 
 from specreduce.calibration_data import load_pypeit_calibration_lines
 from specutils import Spectrum
 from specreduce.extract import BoxcarExtract
-from specreduce.line_matching import match_lines_wcs, find_arc_lines
+from specreduce.line_matching import match_lines_wcs, find_arc_lines, _estimate_baseline
 from specreduce.tracing import FlatTrace
 from specreduce.utils.synth_data import SynthImage
 
@@ -116,3 +117,80 @@ def test_match_lines_wcs(mk_test_data):
         tolerance=5
     )
     assert len(matched_lines) > 1
+
+
+def _mk_pedestal_spectrum(
+    pedestal, line_centers, amplitude=200.0, fwhm=4.0, npix=1000, seed=1
+):
+    """Synthetic 1D arc spectrum: Gaussian lines on top of a (possibly varying) pedestal."""
+    rng = np.random.default_rng(seed)
+    x = np.arange(npix)
+    sigma = fwhm * 0.4247
+    flux = np.broadcast_to(np.asarray(pedestal, dtype=float), (npix,)).copy()
+    for c in line_centers:
+        flux += amplitude * np.exp(-0.5 * ((x - c) / sigma) ** 2)
+    flux += rng.normal(0.0, 1.0, npix)
+    return Spectrum(
+        flux * u.DN,
+        spectral_axis=x * u.pix,
+        uncertainty=StdDevUncertainty(np.ones(npix)),
+    )
+
+
+def test_estimate_baseline_global():
+    """A global sigma-clipped median should recover a flat pedestal under bright lines."""
+    spec = _mk_pedestal_spectrum(50.0, [100, 300, 500, 700, 900])
+    baseline = _estimate_baseline(spec.flux.value)
+    assert baseline.shape == spec.flux.shape
+    np.testing.assert_allclose(baseline, 50.0, atol=1.0)
+
+
+def test_estimate_baseline_windowed():
+    """A windowed estimate should follow a slowly varying pedestal."""
+    npix = 1000
+    ramp = np.linspace(20.0, 120.0, npix)
+    spec = _mk_pedestal_spectrum(ramp, [100, 300, 500, 700, 900], npix=npix)
+    baseline = _estimate_baseline(spec.flux.value, window=100)
+    # Away from the edges (where the interpolation is constant) the estimate tracks the ramp.
+    np.testing.assert_allclose(baseline[100:-100], ramp[100:-100], atol=3.0)
+    # The global estimate cannot follow the ramp.
+    assert np.abs(_estimate_baseline(spec.flux.value) - ramp).max() > 30.0
+
+
+@pytest.mark.parametrize("window", [0, -5, 2.5])
+def test_estimate_baseline_invalid_window(window):
+    with pytest.raises(ValueError, match="window"):
+        _estimate_baseline(np.ones(100), window=window)
+
+
+def test_find_arc_lines_subtract_baseline():
+    """A pedestal above the detection threshold hides every line unless it is subtracted."""
+    centers = [100, 250, 400, 550, 700, 850]
+    spec = _mk_pedestal_spectrum(50.0, centers)
+
+    # Threshold is noise_factor * sigma = 5 < 50, so the whole spectrum is one "line"
+    # and specutils warns that the spectrum is not below the threshold.
+    with pytest.warns(AstropyUserWarning, match="not below the threshold"):
+        lines = find_arc_lines(spec, fwhm=4, noise_factor=5, subtract_baseline=False)
+    assert len(lines) <= 1
+
+    lines = find_arc_lines(spec, fwhm=4, noise_factor=5, subtract_baseline=True)
+    assert len(lines) == len(centers)
+    np.testing.assert_allclose(np.sort(lines["centroid"].value), centers, atol=0.1)
+    # Amplitudes are measured above the baseline, not from zero.
+    np.testing.assert_allclose(lines["amplitude"].value, 200.0, rtol=0.05)
+    # The input spectrum is left untouched.
+    assert np.median(spec.flux.value) > 40.0
+
+
+def test_find_arc_lines_subtract_baseline_windowed():
+    """A ramped pedestal needs the windowed baseline to recover every line."""
+    centers = [100, 250, 400, 550, 700, 850]
+    ramp = np.linspace(20.0, 120.0, 1000)
+    spec = _mk_pedestal_spectrum(ramp, centers)
+
+    lines = find_arc_lines(
+        spec, fwhm=4, noise_factor=5, subtract_baseline=True, baseline_window=100
+    )
+    assert len(lines) == len(centers)
+    np.testing.assert_allclose(np.sort(lines["centroid"].value), centers, atol=0.1)