From 4d58c9fd127a56bec45035d267f63f092c267fe3 Mon Sep 17 00:00:00 2001
From: Hannu Parviainen <hpparvi@gmail.com>
Date: Sat, 5 Sep 2026 12:46:16 +0100
Subject: [PATCH 1/4] - Fixed TiltSolution.resample dropping the input
 uncertainty, mask, and metadata. The uncertainty is now propagated through
 the resampling and returned in the same uncertainty class as the input, the
 mask flags every bin that received a contribution from a masked pixel, and
 the metadata is copied to the resampled NDData.

---
 CHANGES.rst                              |   8 ++
 docs/tilt_correction/tilt_correction.rst |   5 +-
 specreduce/tests/test_tilt_solution.py   | 142 ++++++++++++++++++++++-
 specreduce/tilt_solution.py              |  60 ++++++++--
 4 files changed, 201 insertions(+), 14 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 0464fbdc..656d0e42 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -38,6 +38,14 @@ API Changes
   on instantiation, and its removal has been rescheduled from v2.0 to v1.11.
   Use ``specreduce.wavecal1d.WavelengthCalibration1D`` instead. [#316]
 
+Bug Fixes
+^^^^^^^^^
+
+- ``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
 ^^^^^^^^^^^^^
 
diff --git a/docs/tilt_correction/tilt_correction.rst b/docs/tilt_correction/tilt_correction.rst
index 6398243f..2672f36b 100644
--- a/docs/tilt_correction/tilt_correction.rst
+++ b/docs/tilt_correction/tilt_correction.rst
@@ -262,4 +262,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.
diff --git a/specreduce/tests/test_tilt_solution.py b/specreduce/tests/test_tilt_solution.py
index c59edb62..53ee35fc 100644
--- a/specreduce/tests/test_tilt_solution.py
+++ b/specreduce/tests/test_tilt_solution.py
@@ -1,10 +1,26 @@
+import astropy.units as u
 import numpy as np
 import pytest
+from astropy.io import fits
 from astropy.modeling import models
+from astropy.modeling.models import Shift, Polynomial2D
+from astropy.nddata import (
+    NDData,
+    CCDData,
+    StdDevUncertainty,
+    VarianceUncertainty,
+    InverseVariance,
+)
 
 from specreduce.tilt_solution import TiltSolution, diff_poly2d_x
 
 
+def _linear_ts(ny, nx, shift=0.0, disp_axis=1):
+    """A fit-free tilt solution: detector x = rectified x + shift, for every row."""
+    solution = Shift(0) & Shift(0) | Polynomial2D(1, c0_0=shift, c1_0=1.0)
+    return TiltSolution(solution, disp_axis=disp_axis, image_shape=(ny, nx))
+
+
 def test_diff_poly2d_x_valid_derivative():
     model = models.Polynomial2D(degree=2, c0_0=1, c1_0=2, c2_0=3, c0_1=4, c1_1=5, c0_2=6)
     derivative = diff_poly2d_x(model)
@@ -116,7 +132,124 @@ def test_resample(mk_default_tc, mk_arc_frames):
     tc = mk_default_tc
     tc.find_arc_lines(3.0, 5.0)
     tc.fit(4)
-    tc.solution.resample(arcs[0])
+    result = tc.solution.resample(arcs[0])
+    assert isinstance(result.uncertainty, StdDevUncertainty)
+    assert result.uncertainty.array.shape == result.data.shape
+
+
+@pytest.mark.parametrize(
+    "uncertainty_cls", [StdDevUncertainty, VarianceUncertainty, InverseVariance]
+)
+def test_resample_preserves_uncertainty_type(uncertainty_cls):
+    ny, nx = 4, 12
+    data = np.full((ny, nx), 10.0)
+    image = NDData(data * u.ct, uncertainty=uncertainty_cls(np.full((ny, nx), 4.0)))
+    result = _linear_ts(ny, nx).resample(image)
+    assert isinstance(result.uncertainty, uncertainty_cls)
+    assert result.uncertainty.array.shape == result.data.shape
+
+
+def test_resample_uncertainty_identity():
+    ny, nx = 4, 12
+    data = np.full((ny, nx), 10.0)
+    variance = np.arange(1.0, ny * nx + 1).reshape(ny, nx)
+    image = NDData(data * u.ct, uncertainty=VarianceUncertainty(variance))
+    result = _linear_ts(ny, nx).resample(image)
+    np.testing.assert_allclose(result.data, data)
+    np.testing.assert_allclose(result.uncertainty.array, variance)
+    assert result.uncertainty.unit == u.ct**2
+
+
+def test_resample_uncertainty_half_pixel_shift():
+    """
+    Each rectified bin takes half of two neighboring detector pixels, so the
+    variance is 0.25 + 0.25 of the input variance, not 0.5 + 0.5. The last bin
+    covers half of the last pixel only.
+    """
+    ny, nx = 3, 10
+    data = np.full((ny, nx), 10.0)
+    image = NDData(data * u.ct, uncertainty=StdDevUncertainty(np.full((ny, nx), 3.0)))
+    result = _linear_ts(ny, nx, shift=0.5).resample(image)
+    np.testing.assert_allclose(result.data[:, :-1], 10.0)
+    np.testing.assert_allclose(result.data[:, -1], 5.0)
+    np.testing.assert_allclose(result.uncertainty.array[:, :-1], np.sqrt(4.5), rtol=1e-9)
+    np.testing.assert_allclose(result.uncertainty.array[:, -1], 1.5, rtol=1e-9)
+
+
+@pytest.mark.parametrize(
+    "image",
+    [
+        NDData(np.full((4, 12), 10.0) * u.ct),
+        np.full((4, 12), 10.0),
+        np.full((4, 12), 10.0) * u.ct,
+    ],
+)
+def test_resample_without_uncertainty_returns_none(image):
+    result = _linear_ts(4, 12).resample(image)
+    assert result.uncertainty is None
+
+
+def test_resample_copies_meta():
+    ny, nx = 4, 12
+    data = np.full((ny, nx), 10.0)
+    ts = _linear_ts(ny, nx)
+
+    meta = {"OBJECT": "target", "HISTORY": ["step 1"]}
+    result = ts.resample(NDData(data * u.ct, meta=meta))
+    assert result.meta == meta
+    assert result.meta is not meta
+    result.meta["OBJECT"] = "changed"
+    result.meta["HISTORY"].append("step 2")
+    assert meta == {"OBJECT": "target", "HISTORY": ["step 1"]}
+
+    header = fits.Header([("OBJECT", "target"), ("EXPTIME", 30.0)])
+    result = ts.resample(CCDData(data, unit="ct", meta=header))
+    assert isinstance(result.meta, fits.Header)
+    assert result.meta["OBJECT"] == "target" and result.meta["EXPTIME"] == 30.0
+    result.meta["OBJECT"] = "changed"
+    assert header["OBJECT"] == "target"
+
+    assert len(ts.resample(NDData(data * u.ct)).meta) == 0
+
+
+def test_resample_propagates_mask():
+    ny, nx = 4, 12
+    data = np.full((ny, nx), 10.0)
+    mask = np.zeros((ny, nx), dtype=bool)
+    mask[2, 3] = True
+
+    result = _linear_ts(ny, nx).resample(NDData(data * u.ct, mask=mask), mask_treatment="apply")
+    assert result.mask.dtype == bool
+    np.testing.assert_array_equal(result.mask, mask)
+
+    # a half-pixel shift spreads the masked pixel over the two bins that overlap it
+    result = _linear_ts(ny, nx, shift=0.5).resample(
+        NDData(data * u.ct, mask=mask), mask_treatment="apply"
+    )
+    expected = np.zeros((ny, nx), dtype=bool)
+    expected[2, 2:4] = True
+    np.testing.assert_array_equal(result.mask, expected)
+
+    # fill treatments drop the mask before resampling
+    result = _linear_ts(ny, nx).resample(NDData(data * u.ct, mask=mask), mask_treatment="zero_fill")
+    assert not result.mask.any()
+
+
+def test_resample_disp_axis_0_propagates_arrays():
+    n = 8
+    data = np.arange(1.0, n * n + 1).reshape(n, n)
+    variance = 2.0 * data
+    mask = np.zeros((n, n), dtype=bool)
+    mask[1, 5] = True
+    image = NDData(data * u.ct, uncertainty=VarianceUncertainty(variance), mask=mask)
+
+    result = _linear_ts(n, n, disp_axis=0).resample(image)
+    np.testing.assert_allclose(result.data, data.T)
+    np.testing.assert_allclose(result.uncertainty.array, variance.T)
+    np.testing.assert_array_equal(result.mask, mask.T)
+
+    result = _linear_ts(n, n, disp_axis=0).resample(image, nbins=2 * n)
+    assert result.data.shape == result.uncertainty.array.shape == result.mask.shape == (2 * n, n)
 
 
 @pytest.mark.remote_data
@@ -195,14 +328,15 @@ def test_resample_disp_axis_0(mk_default_tc, mk_arc_frames):
     tc.fit(4)
 
     # Use a square crop so _parse_image works with disp_axis=0
-    from astropy.nddata import NDData
-    import astropy.units as u
     ny = arcs[0].data.shape[0]
-    square = NDData(arcs[0].data[:, :ny] * u.ct)
+    square = NDData(
+        arcs[0].data[:, :ny] * u.ct, uncertainty=StdDevUncertainty(np.full((ny, ny), 5.0))
+    )
 
     ts = tc.solution
     ts.disp_axis = 0
     result = ts.resample(square, nbins=ny)
+    assert result.uncertainty.array.shape == result.data.shape
     # With disp_axis=0, output should be transposed
     assert result.data.shape[0] == ny  # nbins along axis 0
     assert result.data.shape[1] == ny  # cdisp along axis 1
diff --git a/specreduce/tilt_solution.py b/specreduce/tilt_solution.py
index f146bb8f..c056d04d 100644
--- a/specreduce/tilt_solution.py
+++ b/specreduce/tilt_solution.py
@@ -1,4 +1,5 @@
 import warnings
+from copy import deepcopy
 from functools import cached_property
 from typing import Sequence, Literal
 
@@ -7,7 +8,7 @@
 import numpy as np
 from astropy.modeling import models, fitting, Model
 from astropy.modeling.models import Identity, Mapping, Shift, Polynomial2D
-from astropy.nddata import NDData
+from astropy.nddata import NDData, VarianceUncertainty
 from astropy.utils.exceptions import AstropyUserWarning
 from gwcs import coordinate_frames
 from numpy import ndarray
@@ -296,14 +297,34 @@ def resample(
         Returns
         -------
         NDData
-            NDData instance containing the flux values resampled into the uniform grid
-            defined by ``nbins``, ``bounds``, or ``bin_edges``.
+            The flux resampled into the uniform grid defined by ``nbins``, ``bounds``, or
+            ``bin_edges``. If the input carries an uncertainty, the resampled uncertainty is
+            attached in the same uncertainty class as the input. The mask marks every output
+            bin that received a contribution from a masked input pixel, and the metadata is a
+            copy of the input metadata. The WCS is not propagated.
+
+        Notes
+        -----
+        Each output bin is a linear combination of detector pixels,
+        ``F = n * sum_j k_j f_j`` with ``k_j`` the fractional pixel overlap times the
+        Jacobian of the transformation and ``n`` the per-row flux-conservation factor.
+        The variance is propagated as ``Var = n**2 * sum_j k_j**2 var_j``, assuming
+        independent pixel noise and treating ``n`` and the Jacobian as deterministic.
         """
 
+        # The metadata and the presence of an uncertainty are read from the input itself:
+        # parse_image drops the metadata and fabricates a unit variance for bare arrays.
+        meta = deepcopy(getattr(flux, "meta", None))
+        has_uncertainty = getattr(flux, "uncertainty", None) is not None
+
         # TODO: In the future, we want to make sure that we don't copy the data unless absolutely
         # necessary.
         im = parse_image(flux, disp_axis=self.disp_axis, mask_treatment=mask_treatment)
         flux = im.flux.value
+        mask = im.mask.astype(bool)
+        if has_uncertainty:
+            uncertainty_type = type(im.uncertainty)
+            variance = im.uncertainty.represent_as(VarianceUncertainty).array
 
         ny, nx = flux.data.shape
         ypix = np.arange(ny)
@@ -316,6 +337,8 @@ def resample(
         bin_edge_w = bin_edges_det - bin_edge_ix
 
         resampled_flux = np.zeros((ny, nbins))
+        resampled_variance = np.zeros((ny, nbins)) if has_uncertainty else None
+        resampled_mask = np.zeros((ny, nbins), dtype=bool)
         weights = np.zeros((ny, nx))
 
         # Calculate the derivative of the tilt-corrected space -> detector space transformation with
@@ -343,9 +366,11 @@ def resample(
             # the tilt-corrected flux is the detector flux in that pixel, scaled by the width of the
             # tilt-corrected bin in detector coordinates and the derivative dtdx.
             if m.any():
-                resampled_flux[:, i] = (
-                    (bin_edges_det[:, i + 1] - bin_edges_det[:, i]) * flux[ys, i1] * dtdx[ys, i1]
-                )
+                k = (bin_edges_det[:, i + 1] - bin_edges_det[:, i]) * dtdx[ys, i1]
+                resampled_flux[:, i] = k * flux[ys, i1]
+                resampled_mask[:, i] = mask[ys, i1]
+                if has_uncertainty:
+                    resampled_variance[:, i] = k**2 * variance[ys, i1]
 
             # For rows where the tilt-corrected bin spans multiple detector pixels, calculate the
             # tilt-corrected flux as a weighted sum of the detector flux, multiplied by dtdx,
@@ -358,11 +383,28 @@ def resample(
                 w[(ixc > i1[:, None]) & (ixc < i2[:, None])] = 1
                 w[ys, i1 - imin] = 1.0 - bin_edge_w[:, i]
                 w[ys, i2 - imin] = bin_edge_w[:, i + 1]
-                resampled_flux[~m, i] = (flux[~m, imin:imax] * dtdx[~m, imin:imax] * w[~m]).sum(1)
+                k = dtdx[~m, imin:imax] * w[~m]
+                resampled_flux[~m, i] = (flux[~m, imin:imax] * k).sum(1)
+                resampled_mask[~m, i] = (mask[~m, imin:imax] & (w[~m] > 0)).any(1)
+                if has_uncertainty:
+                    resampled_variance[~m, i] = (variance[~m, imin:imax] * k**2).sum(1)
 
         # Apply the normalization factor to conserve flux
         resampled_flux *= n[:, None]
+        if has_uncertainty:
+            resampled_variance *= n[:, None] ** 2
+
         if self.disp_axis == 0:
             resampled_flux = resampled_flux.T
-
-        return NDData(resampled_flux * im.unit)
+            resampled_mask = resampled_mask.T
+            if has_uncertainty:
+                resampled_variance = resampled_variance.T
+
+        uncertainty = None
+        if has_uncertainty:
+            uncertainty = VarianceUncertainty(resampled_variance * im.unit**2).represent_as(
+                uncertainty_type
+            )
+        return NDData(
+            resampled_flux * im.unit, uncertainty=uncertainty, mask=resampled_mask, meta=meta
+        )

From 778ae5bd6d4d8b4086042040dc43da3ac9562a4b Mon Sep 17 00:00:00 2001
From: Hannu Parviainen <hpparvi@gmail.com>
Date: Tue, 8 Sep 2026 16:45:04 +0100
Subject: [PATCH 2/4] Added baseline subtraction to arc line detection.
 `line_matching.find_arc_lines` gained `subtract_baseline` and
 `baseline_window` options that estimate the baseline flux with a
 sigma-clipped median (globally, or in running windows with linear edge
 extrapolation) and remove it before the threshold-based line detection.
 `TiltCorrection.find_arc_lines` exposes the same options with subtraction on
 by default, so it no longer requires background-subtracted arc frames to
 detect any lines.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_012Pfq3Gc7Boe9TKDs7DVE64
---
 CHANGES.rst                              |   7 ++
 docs/tilt_correction/tilt_correction.rst |  13 ++-
 specreduce/line_matching.py              | 108 ++++++++++++++++++++++-
 specreduce/tests/test_line_matching.py   |  80 ++++++++++++++++-
 specreduce/tests/test_tilt_correction.py |  32 +++++++
 specreduce/tilt_correction.py            |  38 +++++++-
 6 files changed, 269 insertions(+), 9 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 656d0e42..5c95aa79 100644
--- a/CHANGES.rst
+++ b/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
 ^^^^^^^^^^^
 
diff --git a/docs/tilt_correction/tilt_correction.rst b/docs/tilt_correction/tilt_correction.rst
index 2672f36b..0a7b5108 100644
--- a/docs/tilt_correction/tilt_correction.rst
+++ b/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
 *************************
diff --git a/specreduce/line_matching.py b/specreduce/line_matching.py
index f8f1b807..34bcdc55 100644
--- a/specreduce/line_matching.py
+++ b/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 = False,
+    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 = []
diff --git a/specreduce/tests/test_line_matching.py b/specreduce/tests/test_line_matching.py
index a96d0675..db087977 100644
--- a/specreduce/tests/test_line_matching.py
+++ b/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)
diff --git a/specreduce/tests/test_tilt_correction.py b/specreduce/tests/test_tilt_correction.py
index 71010004..fd0a0fd2 100644
--- a/specreduce/tests/test_tilt_correction.py
+++ b/specreduce/tests/test_tilt_correction.py
@@ -1,3 +1,5 @@
+from copy import deepcopy
+
 import matplotlib.pyplot as plt
 import numpy as np
 import pytest
@@ -105,3 +107,33 @@ def test_plot_fit_quality_with_rlim(mk_default_tc):
     tc.fit(4)
     tc.plot_fit_quality(rlim=(-1, 1))
     plt.close("all")
+
+
+@pytest.mark.remote_data
+def test_find_lines_bright_background(mk_arc_frames):
+    """A pedestal above the detection threshold hides the lines unless it is subtracted."""
+    arcs = mk_arc_frames
+    reference = TiltCorrection(arc_frames=arcs, cdisp_ref_pixel=64, n_cdisp_samples=8)
+    reference.find_arc_lines(3.0, 5.0)
+
+    # Add a pedestal of 40 sigma (uncertainty is 5) on top of the arc lines.
+    bright = [deepcopy(arc) for arc in arcs]
+    for arc in bright:
+        arc.data = arc.data + 200.0
+    tc = TiltCorrection(arc_frames=bright, cdisp_ref_pixel=64, n_cdisp_samples=8)
+
+    # Without baseline subtraction the whole row is above threshold: one "line" at most.
+    tc.find_arc_lines(3.0, 5.0, subtract_baseline=False)
+    assert all(lines.size <= 1 for lines in tc._lines_ref)
+
+    # With the default baseline subtraction the result is the same as for the faint frames.
+    tc.find_arc_lines(3.0, 5.0)
+    for ref_lines, lines, det_x in zip(reference._lines_ref, tc._lines_ref, tc._samples_det_x):
+        assert lines.size > 1
+        np.testing.assert_allclose(lines, ref_lines, atol=0.05)
+        assert det_x.size >= 0.8 * lines.size * tc.cd_samples.size
+
+    # A running baseline estimate works too.
+    tc.find_arc_lines(3.0, 5.0, baseline_window=64)
+    for ref_lines, lines in zip(reference._lines_ref, tc._lines_ref):
+        np.testing.assert_allclose(lines, ref_lines, atol=0.2)
diff --git a/specreduce/tilt_correction.py b/specreduce/tilt_correction.py
index 4f6da96d..c62d2302 100644
--- a/specreduce/tilt_correction.py
+++ b/specreduce/tilt_correction.py
@@ -141,13 +141,25 @@ def __init__(
 
         self.solution: TiltSolution | None = None
 
-    def find_arc_lines(self, fwhm: float, noise_factor: float = 5.0) -> None:
+    def find_arc_lines(
+        self,
+        fwhm: float,
+        noise_factor: float = 5.0,
+        subtract_baseline: bool = True,
+        baseline_window: int | None = None,
+    ) -> None:
         """Find arc lines from the provided arc frames for all cross-dispersion samples.
 
         This method locates spectral arc lines from the provided arc frames, calculates
         their centroids, and organizes them into reference lists and sample arrays
         for further analysis.
 
+        The arc frames are used as they are, so by default the baseline (background)
+        flux of each row is estimated with a sigma-clipped median and removed before
+        the line detection. The detection thresholds the flux against
+        ``noise_factor × uncertainty``, so a pedestal above that level would otherwise
+        swallow all the lines.
+
         Parameters
         ----------
         fwhm
@@ -155,6 +167,14 @@ def find_arc_lines(self, fwhm: float, noise_factor: float = 5.0) -> None:
             by the line-finding algorithm.
         noise_factor
             A multiplier for noise thresholding in the line-finding process.
+        subtract_baseline
+            Estimate and subtract the baseline flux of each row before the line
+            detection. Set to ``False`` for arc frames that are already
+            background-subtracted.
+        baseline_window
+            Width in pixels of the chunks used for a running baseline estimate that
+            can follow a slowly varying background. If ``None``, a single global
+            median is subtracted from each row.
         """
         self._arc_spectra = []
         self._samples_rec_x = []
@@ -174,7 +194,13 @@ def find_arc_lines(self, fwhm: float, noise_factor: float = 5.0) -> None:
                     d.data[self.ref_pixel[0]] * d.unit,
                     uncertainty=d.uncertainty[self.ref_pixel[0]].represent_as(StdDevUncertainty),
                 )
-                lines = find_arc_lines(spectrum, fwhm, noise_factor=noise_factor)
+                lines = find_arc_lines(
+                    spectrum,
+                    fwhm,
+                    noise_factor=noise_factor,
+                    subtract_baseline=subtract_baseline,
+                    baseline_window=baseline_window,
+                )
                 self._lines_ref.append(lines["centroid"].value)
 
                 # Find the line centroids for the sample rows
@@ -183,7 +209,13 @@ def find_arc_lines(self, fwhm: float, noise_factor: float = 5.0) -> None:
                         d.data[s] * d.unit,
                         uncertainty=d.uncertainty[s].represent_as(StdDevUncertainty),
                     )
-                    lines = find_arc_lines(spectrum, fwhm, noise_factor=noise_factor)
+                    lines = find_arc_lines(
+                        spectrum,
+                        fwhm,
+                        noise_factor=noise_factor,
+                        subtract_baseline=subtract_baseline,
+                        baseline_window=baseline_window,
+                    )
                     samples_x[i].append(lines["centroid"].value)
                     samples_y[i].append(np.full(len(lines), s))
                     self._arc_spectra[i].append(spectrum)

From f3c87e2cb450a4e42440412dc772fff17333650b Mon Sep 17 00:00:00 2001
From: Hannu Parviainen <hpparvi@gmail.com>
Date: Tue, 8 Sep 2026 17:10:53 +0100
Subject: [PATCH 3/4] - Improved `TiltCorrection.find_arc_lines` stability.

---
 CHANGES.rst                              |  7 ++++
 specreduce/tests/test_tilt_correction.py | 51 ++++++++++++++++++++++++
 specreduce/tilt_correction.py            | 25 ++++++++++++
 3 files changed, 83 insertions(+)

diff --git a/CHANGES.rst b/CHANGES.rst
index 5c95aa79..136683c1 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -48,6 +48,13 @@ API Changes
 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``.
diff --git a/specreduce/tests/test_tilt_correction.py b/specreduce/tests/test_tilt_correction.py
index fd0a0fd2..3b6a4d83 100644
--- a/specreduce/tests/test_tilt_correction.py
+++ b/specreduce/tests/test_tilt_correction.py
@@ -137,3 +137,54 @@ def test_find_lines_bright_background(mk_arc_frames):
     tc.find_arc_lines(3.0, 5.0, baseline_window=64)
     for ref_lines, lines in zip(reference._lines_ref, tc._lines_ref):
         np.testing.assert_allclose(lines, ref_lines, atol=0.2)
+
+
+@pytest.mark.remote_data
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        dict(cdisp_ref_pixel=64, cdisp_samples=[10, 60, 200]),
+        dict(cdisp_ref_pixel=64, cdisp_samples=[-5, 60, 100]),
+        dict(cdisp_ref_pixel=64, cdisp_sample_lims=(0, 300)),
+        dict(cdisp_ref_pixel=300, n_cdisp_samples=4),
+    ],
+)
+def test_find_lines_rejects_samples_outside_frame(mk_arc_frames, kwargs):
+    """The reference row and all cross-dispersion samples must lie inside the image frame."""
+    tc = TiltCorrection(mk_arc_frames, **kwargs)
+    with pytest.raises(ValueError, match="outside the image frame"):
+        tc.find_arc_lines(3.0, 5.0)
+
+
+@pytest.mark.remote_data
+def test_find_lines_drops_centroids_outside_frame(mk_default_tc, monkeypatch):
+    """Fitted line centroids that fall off the frame along the dispersion axis are dropped."""
+    import specreduce.tilt_correction as tcmod
+    from astropy.table import vstack
+
+    tc = mk_default_tc
+    nx = tc.arc_frames[0].data.shape[1]
+
+    tc.find_arc_lines(3.0, 5.0)
+    expected_ref = [lines.copy() for lines in tc._lines_ref]
+    expected_det = [x.copy() for x in tc._samples_det_x]
+
+    # Make every line search also return three bogus lines: one off each end of the
+    # frame and one with a failed (NaN) fit.
+    real_find_arc_lines = tcmod.find_arc_lines
+
+    def bogus_find_arc_lines(*args, **kwargs):
+        lines = real_find_arc_lines(*args, **kwargs)
+        extra = lines[:3].copy()
+        extra["centroid"] = [-3.0, nx + 2.0, np.nan] * extra["centroid"].unit
+        return vstack([lines, extra])
+
+    monkeypatch.setattr(tcmod, "find_arc_lines", bogus_find_arc_lines)
+    tc.find_arc_lines(3.0, 5.0)
+
+    for got, want in zip(tc._lines_ref, expected_ref):
+        np.testing.assert_array_equal(got, want)
+    for got, want in zip(tc._samples_det_x, expected_det):
+        np.testing.assert_array_equal(got, want)
+    for x in tc._samples_rec_x + tc._samples_det_x:
+        assert np.all((x >= 0) & (x < nx))
diff --git a/specreduce/tilt_correction.py b/specreduce/tilt_correction.py
index c62d2302..015f97e9 100644
--- a/specreduce/tilt_correction.py
+++ b/specreduce/tilt_correction.py
@@ -5,6 +5,7 @@
 import numpy as np
 from astropy.modeling import models, fitting
 from astropy.nddata import StdDevUncertainty, NDData
+from astropy.table import QTable
 from numpy import ndarray, repeat, tile
 from scipy.optimize import minimize
 from scipy.spatial import KDTree
@@ -141,6 +142,11 @@ def __init__(
 
         self.solution: TiltSolution | None = None
 
+    def _lines_inside_frame(self, lines: QTable) -> QTable:
+        """Drop detected lines whose centroid is not finite or falls off the frame."""
+        x = lines["centroid"].value
+        return lines[np.isfinite(x) & (x >= 0) & (x < self._nx)]
+
     def find_arc_lines(
         self,
         fwhm: float,
@@ -175,7 +181,24 @@ def find_arc_lines(
             Width in pixels of the chunks used for a running baseline estimate that
             can follow a slowly varying background. If ``None``, a single global
             median is subtracted from each row.
+
+        Detected lines whose fitted centroid falls outside the image frame along the
+        dispersion axis (or is not finite) are discarded, as they come from failed fits.
+
+        Raises
+        ------
+        ValueError
+            If the reference row or any cross-dispersion sample lies outside the
+            image frame.
         """
+        rows = np.append(self.cd_samples, self.ref_pixel[0])
+        outside = rows[(rows < 0) | (rows >= self._ny)]
+        if outside.size > 0:
+            raise ValueError(
+                f"Cross-dispersion samples {outside.tolist()} lie outside the image frame; "
+                f"samples must satisfy 0 <= sample < {self._ny}."
+            )
+
         self._arc_spectra = []
         self._samples_rec_x = []
         self._lines_ref = []
@@ -201,6 +224,7 @@ def find_arc_lines(
                     subtract_baseline=subtract_baseline,
                     baseline_window=baseline_window,
                 )
+                lines = self._lines_inside_frame(lines)
                 self._lines_ref.append(lines["centroid"].value)
 
                 # Find the line centroids for the sample rows
@@ -216,6 +240,7 @@ def find_arc_lines(
                         subtract_baseline=subtract_baseline,
                         baseline_window=baseline_window,
                     )
+                    lines = self._lines_inside_frame(lines)
                     samples_x[i].append(lines["centroid"].value)
                     samples_y[i].append(np.full(len(lines), s))
                     self._arc_spectra[i].append(spectrum)

From 41f85a512e2a6370438f25099d1a8efb96bcb9cd Mon Sep 17 00:00:00 2001
From: Hannu Parviainen <hpparvi@gmail.com>
Date: Thu, 10 Sep 2026 11:34:56 +0100
Subject: [PATCH 4/4] - Changed `find_arc_lines` to set `subtract_baseline`
 True by default.

---
 specreduce/line_matching.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/specreduce/line_matching.py b/specreduce/line_matching.py
index 34bcdc55..fd755af3 100644
--- a/specreduce/line_matching.py
+++ b/specreduce/line_matching.py
@@ -95,7 +95,7 @@ def find_arc_lines(
     fwhm: float | u.Quantity = 5.0 * u.pix,
     window: float = 3.0,
     noise_factor: float = 5.0,
-    subtract_baseline: bool = False,
+    subtract_baseline: bool = True,
     baseline_window: int | None = None,
 ) -> QTable:
     """