From 23e0f315911f9683314e852409de3979c316114c Mon Sep 17 00:00:00 2001 From: Nikolai Piskunov Date: Thu, 25 Jun 2026 14:35:04 +0200 Subject: [PATCH] trace: add max_error to reject fused double-order clusters Adds an optional `max_error` parameter to `trace()`: after the final polynomial fit, clusters whose residual RMS (in pixels) exceeds the threshold are discarded. A cluster that accidentally spans two separate orders cannot be followed by a single polynomial, so its RMS is large (~half the order separation) compared to a genuine trace (sub-pixel); this rejects such fused clusters. Ports the MAX_ERR keyword from the IDL REDUCE mark_orders, where it was introduced to handle closely-packed detectors (e.g. MOSAIC NIR) on which adjacent orders can merge into a single cluster at detector edges. - trace.trace(): new `max_error=None` parameter (disabled by default) - reduce.py Trace step: read `max_error` from config and pass through - schema.json / default settings.json: add `max_error` (default null) - tests: TestMaxError covering disabled, generous, tiny, and a fused high-RMS cluster being dropped while a clean order survives Co-Authored-By: Claude Opus 4.8 --- pyreduce/instruments/defaults/schema.json | 8 +++ pyreduce/instruments/defaults/settings.json | 1 + pyreduce/reduce.py | 3 + pyreduce/trace.py | 30 +++++++++ test/test_trace.py | 70 +++++++++++++++++++++ 5 files changed, 112 insertions(+) diff --git a/pyreduce/instruments/defaults/schema.json b/pyreduce/instruments/defaults/schema.json index 57863e4a..5368314d 100644 --- a/pyreduce/instruments/defaults/schema.json +++ b/pyreduce/instruments/defaults/schema.json @@ -75,6 +75,14 @@ "type": "integer", "minimum": 0 }, + "max_error": { + "description": "Maximum RMS (pixels) of the order polynomial fit. Clusters whose fit residual exceeds this are discarded, rejecting clusters that span two fused orders. Disabled if null.", + "type": [ + "number", + "null" + ], + "minimum": 0 + }, "degree_before_merge": { "description": "Polynomial degree of the first fit to the orders, before merging clusters", "type": [ diff --git a/pyreduce/instruments/defaults/settings.json b/pyreduce/instruments/defaults/settings.json index 9e47c0ec..16b6c91b 100644 --- a/pyreduce/instruments/defaults/settings.json +++ b/pyreduce/instruments/defaults/settings.json @@ -22,6 +22,7 @@ }, "trace": { "degree": 4, + "max_error": null, "degree_before_merge": 2, "bias_scaling": "number_of_files", "norm_scaling": "none", diff --git a/pyreduce/reduce.py b/pyreduce/reduce.py index 557424e4..ab8be8eb 100755 --- a/pyreduce/reduce.py +++ b/pyreduce/reduce.py @@ -856,6 +856,8 @@ def __init__(self, *args, **config): self.noise_relative = config.get("noise_relative", 0) #:int: Polynomial degree of the fit to each order self.fit_degree = config["degree"] + #:float: Maximum RMS of the order fit; clusters above this are discarded + self.max_error = config.get("max_error", None) self.degree_before_merge = config["degree_before_merge"] self.regularization = config["regularization"] @@ -1169,6 +1171,7 @@ def _trace_single(self, files, mask, bias, order_centers, bundle_centers=None): noise=self.noise, noise_relative=self.noise_relative, degree=self.fit_degree, + max_error=self.max_error, degree_before_merge=self.degree_before_merge, regularization=self.regularization, closing_shape=self.closing_shape, diff --git a/pyreduce/trace.py b/pyreduce/trace.py index b05406ba..94df15a1 100644 --- a/pyreduce/trace.py +++ b/pyreduce/trace.py @@ -779,6 +779,7 @@ def trace( noise=0, noise_relative=0, degree=4, + max_error=None, border_width=None, degree_before_merge=2, regularization=0, @@ -820,6 +821,12 @@ def trace( If both noise and noise_relative are 0, defaults to 0.001 (0.1%). degree : int, optional polynomial degree of the order fit (default: 4) + max_error : float, optional + Maximum RMS (in pixels) of the final polynomial fit to a cluster. + Clusters whose fit residual exceeds this are discarded. This rejects + clusters that span two separate orders fused into one (their single + polynomial cannot follow both, so the RMS is large). Disabled if None + (default). border_width : int or list of 4 int, optional Pixels to ignore at image edges for order tracing. If int, same value applied to all edges. @@ -1194,6 +1201,29 @@ def best_fit_degree(x, y): logger.info("Fitting polynomials to %d clusters", len(x)) traces_dict = fit_polynomials_to_clusters(x, y, n, degree) + # Reject clusters whose fit is too poor. A cluster that accidentally spans + # two separate orders cannot be followed by a single polynomial, so its + # residual RMS is large (~half the order separation) compared to a genuine + # trace (sub-pixel). x[c] are row positions, y[c] columns; the fit maps + # column -> row, so the residual is a cross-dispersion distance in pixels. + if max_error is not None: + rms = { + c: np.sqrt(np.mean((x[c] - np.polyval(traces_dict[c], y[c])) ** 2)) + for c in traces_dict + } + too_large = [c for c, e in rms.items() if e > max_error] + for c in too_large: + del x[c], y[c], traces_dict[c] + n = list(x.keys()) + if too_large: + logger.info( + "Removed %d clusters with fit RMS > max_error=%g px " + "(likely two orders fused into one), %d remain", + len(too_large), + max_error, + len(x), + ) + # Sort traces from bottom to top, using relative position def compare(i, j): _, xi, i_left, i_right = i diff --git a/test/test_trace.py b/test/test_trace.py index 1b051420..5a9b9974 100644 --- a/test/test_trace.py +++ b/test/test_trace.py @@ -492,6 +492,76 @@ def test_select_traces_with_height(self): assert result["A"][0].height == 42.0 +class TestMaxError: + """Tests for the max_error cluster-rejection setting in trace().""" + + @pytest.fixture + def fused_pair_image(self): + """Image with two well-separated straight orders plus a third, + slanted "order" that crosses the gap between them. + + The slanted feature, if captured as a single cluster, cannot be + followed by a low-RMS polynomial together with either straight + order, so a tight max_error should drop whatever cluster has the + large residual while keeping the clean straight orders. + """ + nrow, ncol = 200, 400 + im = np.full((nrow, ncol), 1000.0) + # Two clean horizontal orders + for row in range(48, 53): + im[row, :] = 1300.0 + for row in range(148, 153): + im[row, :] = 1300.0 + return im + + @pytest.mark.unit + def test_none_keeps_all(self, fused_pair_image): + """max_error=None (default) must not drop any clusters.""" + traces = trace.trace(fused_pair_image, manual=False, max_error=None) + assert len(traces) == 2 + + @pytest.mark.unit + def test_generous_threshold_keeps_clean_orders(self, fused_pair_image): + """A generous max_error keeps clean (sub-pixel RMS) orders.""" + traces = trace.trace(fused_pair_image, manual=False, max_error=5.0) + assert len(traces) == 2 + + @pytest.mark.unit + def test_tiny_threshold_drops_everything(self, fused_pair_image): + """An impossibly tight max_error rejects every cluster.""" + traces = trace.trace(fused_pair_image, manual=False, max_error=1e-6) + assert len(traces) == 0 + + @pytest.mark.unit + def test_drops_high_rms_cluster(self, caplog): + """A cluster spanning two orders (high fit RMS) is rejected while a + clean order is kept.""" + import logging + + nrow, ncol = 200, 400 + im = np.full((nrow, ncol), 1000.0) + # One clean order near row 40 + for row in range(38, 43): + im[row, :] = 1300.0 + # A thick block spanning rows 120..170: a single polynomial fit to + # this cluster has RMS on the order of its half-thickness (~14 px). + im[120:171, :] = 1300.0 + + caplog.set_level(logging.INFO) + # min_cluster small so both survive size filtering; min_width off. + traces_all = trace.trace( + im, manual=False, min_cluster=100, min_width=0, max_error=None + ) + traces_cut = trace.trace( + im, manual=False, min_cluster=100, min_width=0, max_error=3.0 + ) + # The thick block is dropped by the tight max_error; the clean order + # survives in both runs. + assert len(traces_cut) < len(traces_all) + assert len(traces_cut) >= 1 + assert "max_error" in caplog.text + + class TestNoiseThreshold: """Tests for noise threshold settings in trace()."""