Expose rolling_ball's hidden interpolation orders as proper kwargs

decompose_rolling_ball inlined a down→process→up round-trip with order=1
on the way down and order=3 on the way up, both hardcoded. That asymmetry
shaped every radius>150 result in the existing 39,731-combination sweep
but was invisible to any parameter recording — the redundancy analysis
could not tell whether two rolling_ball outputs differed because of radius
or because of the fixed-but-unrecorded interpolation pair.

Add roundtrip_resample(arr, factor, down_order, up_order, transform=None)
in src/utils/preprocessing.py as a reusable helper. Refactor
decompose_rolling_ball to take down_order, up_order, target_effective_radius,
and large_radius_threshold as kwargs, with defaults that preserve current
behavior bit-for-bit so existing sweep results stay valid. Register
down_order and up_order in the decomposition registry so they show up in
any future sweep's recorded params.

All 166 tests pass. Verified directly that the round-trip path (radius=200)
produces the expected shape and that the newly exposed kwargs change the
output (L2 distance ~1126 between default and symmetric orders).

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

Author: bshepp

src/decomposition/methods_extended.py

@@ -460,48 +460,75 @@ def decompose_anisotropic_diffusion(
 @register_decomposition(
     name="rolling_ball",
     category="morphological",
-    default_params={"radius": 50},
-    param_ranges={"radius": [10, 25, 50, 100, 200]},
+    default_params={
+        "radius": 50,
+        "down_order": 1,
+        "up_order": 3,
+        "target_effective_radius": 75,
+        "large_radius_threshold": 150,
+    },
+    param_ranges={
+        "radius": [10, 25, 50, 100, 200],
+        "down_order": [0, 1, 2, 3],
+        "up_order": [1, 3, 5],
+    },
     preserves="features smaller than ball radius",
     destroys="background curvature, large-scale variation",
 )
-def decompose_rolling_ball(dem: np.ndarray, radius: int = 50) -> tuple:
+def decompose_rolling_ball(
+    dem: np.ndarray,
+    radius: int = 50,
+    down_order: int = 1,
+    up_order: int = 3,
+    target_effective_radius: int = 75,
+    large_radius_threshold: int = 150,
+) -> tuple:
     """
     Rolling ball background subtraction.
 
     Simulates a ball rolling under the surface. Common in microscopy.
     Effective for removing large-scale curvature while preserving local features.
 
-    For large radii (>150), uses a downsampling approach to reduce memory usage
-    while maintaining equivalent results for trend extraction.
-    """
-    from scipy.ndimage import zoom
+    For large radii (> large_radius_threshold), runs on a downsampled copy of
+    the DEM to bound memory and CPU. The downsample factor is chosen so the
+    effective ball radius on the small grid is ~target_effective_radius.
+
+    Args:
+        radius: Rolling ball radius in pixels. Smaller = finer features kept.
+        down_order: scipy.ndimage.zoom interpolation order used on the
+            downsample step of the round-trip. Default 1 (linear). Previously
+            hardcoded; exposed so it can be swept and recorded.
+        up_order: scipy.ndimage.zoom interpolation order used on the upsample
+            step of the round-trip. Default 3 (cubic). Asymmetric with
+            down_order by historical default — exposed so callers can pick
+            symmetric pairs if they want.
+        target_effective_radius: When downsampling, the ball radius is rescaled
+            so that the small-grid radius is approximately this value. Lower =
+            coarser intermediate grid = faster but blurrier trend. Default 75.
+        large_radius_threshold: Radii above this trigger the downsample
+            round-trip; radii at or below run directly on the full grid.
+            Default 150.
+    """
+    from ..utils.preprocessing import roundtrip_resample
 
     dem_filled = fill_nans(dem)
     original_shape = dem_filled.shape
 
-    # For large radii, use downsampling to reduce memory
-    # Threshold chosen based on memory constraints
-    large_radius_threshold = 150
-
     if radius > large_radius_threshold:
-        # Downsample factor: aim for effective radius of ~75
-        downsample_factor = radius / 75
-
-        # Downsample the DEM
-        dem_small = zoom(dem_filled, 1.0 / downsample_factor, order=1)
+        downsample_factor = radius / target_effective_radius
         effective_radius = int(radius / downsample_factor)
 
-        # Run rolling ball on smaller DEM
-        trend_small = _rolling_ball_core(dem_small, effective_radius)
-
-        # Upsample the trend back to original size
-        trend = zoom(trend_small, downsample_factor, order=3)
-
-        # Handle size mismatch from zoom rounding
-        if trend.shape != original_shape:
-            # Crop or pad to match
-            trend = _match_shape(trend, original_shape)
+        # Run rolling ball on a coarse copy, then round-trip back to full size.
+        # The down/up interpolation orders are deliberately separable kwargs
+        # so the asymmetry shows up in any parameter sweep that records them.
+        trend = roundtrip_resample(
+            dem_filled,
+            factor=downsample_factor,
+            down_order=down_order,
+            up_order=up_order,
+            transform=lambda small: _rolling_ball_core(small, effective_radius),
+            target_shape=original_shape,
+        )
     else:
         trend = _rolling_ball_core(dem_filled, radius)
 

src/utils/preprocessing.py

@@ -5,6 +5,7 @@
 """
 
 import numpy as np
+from scipy.ndimage import zoom
 
 
 def fill_nans(dem: np.ndarray) -> np.ndarray:
@@ -20,3 +21,56 @@ def fill_nans(dem: np.ndarray) -> np.ndarray:
         Array with NaN values replaced by the mean of non-NaN values.
     """
     return np.nan_to_num(dem, nan=np.nanmean(dem))
+
+
+def roundtrip_resample(
+    arr: np.ndarray,
+    factor: float,
+    down_order: int = 1,
+    up_order: int = 3,
+    transform=None,
+    target_shape: tuple = None,
+) -> np.ndarray:
+    """Downsample arr by 1/factor, optionally transform on the coarse grid,
+    then upsample back by factor.
+
+    Used by decompositions that operate on a coarser grid for memory or
+    speed reasons (e.g. rolling_ball at large radii). Returning the
+    round-tripped array as a single utility makes the two interpolation
+    orders visible to callers — they used to be hardcoded inside callers,
+    invisible to any parameter sweep or reproducibility log.
+
+    Args:
+        arr: 2-D input array.
+        factor: Downsample factor (>1 means coarser intermediate grid).
+        down_order: scipy.ndimage.zoom interpolation order on the way down.
+            Default 1 (linear).
+        up_order: scipy.ndimage.zoom interpolation order on the way up.
+            Default 3 (cubic). down_order and up_order are deliberately
+            independent — callers can pick asymmetric pairs that match
+            their tolerance for smoothing vs. ringing in each direction.
+        transform: Optional callable applied to the downsampled array
+            before upsampling, signature transform(small_arr) -> small_arr.
+            Lets callers do work on the coarse grid (e.g. morphological
+            opening at a smaller effective radius) while keeping the
+            resample bookkeeping in one place.
+        target_shape: Optional (h, w) to crop/pad the output to. Useful
+            when zoom rounding produces an off-by-one mismatch with the
+            input shape and the caller needs an exact size match.
+
+    Returns:
+        Resampled 2-D array. Shape matches target_shape if given,
+        otherwise approximately matches the original (may differ by 1
+        pixel per axis due to zoom rounding at non-integer factors).
+    """
+    small = zoom(arr, 1.0 / factor, order=down_order)
+    if transform is not None:
+        small = transform(small)
+    upscaled = zoom(small, factor, order=up_order)
+    if target_shape is not None and upscaled.shape != target_shape:
+        out = np.zeros(target_shape, dtype=upscaled.dtype)
+        min_h = min(upscaled.shape[0], target_shape[0])
+        min_w = min(upscaled.shape[1], target_shape[1])
+        out[:min_h, :min_w] = upscaled[:min_h, :min_w]
+        return out
+    return upscaled