diff --git a/CHANGES.rst b/CHANGES.rst index 0464fbdc..136683c1 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 ^^^^^^^^^^^ @@ -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 ^^^^^^^^^^^^^ diff --git a/docs/tilt_correction/tilt_correction.rst b/docs/tilt_correction/tilt_correction.rst index 6398243f..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 ************************* @@ -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. diff --git a/specreduce/line_matching.py b/specreduce/line_matching.py index f8f1b807..fd755af3 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 = 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 = [] 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..3b6a4d83 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,84 @@ 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) + + +@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/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_correction.py b/specreduce/tilt_correction.py index 4f6da96d..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,13 +142,30 @@ def __init__( self.solution: TiltSolution | None = None - def find_arc_lines(self, fwhm: float, noise_factor: float = 5.0) -> 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, + 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,7 +173,32 @@ 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. + + 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 = [] @@ -174,7 +217,14 @@ 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, + ) + lines = self._lines_inside_frame(lines) self._lines_ref.append(lines["centroid"].value) # Find the line centroids for the sample rows @@ -183,7 +233,14 @@ 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, + ) + 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) 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 + )