import numpy as np
import time
import imcombiners as imc
import imcombiners.kernels as imck
# Helpers used by the grow timing example below.
def bench_case(label, func, repeats=5):
func() # warm-up
elapsed = []
last = None
for _ in range(repeats):
start = time.perf_counter()
last = func()
elapsed.append(time.perf_counter() - start)
return label, float(np.median(elapsed) * 1_000.0), last
def median_only():
cmb = imc.Combiner(bench)
return cmb.combine("median"), None
def sigclip_median(grow):
c = imc.Combiner(bench).reject(
imc.SigClip(sigma=3.0, maxiters=5, ddof=0, nkeep=0),
grow=grow,
)
return c.combine("median"), c.mask_rej05 Rejection Methods
imcombiners has five rejection algorithms. This tutorial uses Standard Combiner for normal use. It uses Chained Combiner only when an example needs to inspect diagnostics. Direct kernel calls appear at the end.
1. Test Stack
We build a controlled stack: Gaussian noise plus a handful of extreme outliers.
rng = np.random.default_rng(20250311)
n, ny, nx = 25, 8, 8
true_val = 100.0
sigma_noise = 5.0
stack = rng.normal(true_val, sigma_noise, (n, ny, nx)).astype("float32")
# Inject outliers at known locations
stack[3, 2, 2] = 500.0 # bright cosmic ray
stack[7, 5, 4] = 480.0 # another cosmic ray
stack[14, 0, 6] = -200.0 # negative spike
stack[0] += 40.0 # one systematically hot frameFor production reductions, pass a Rejector object to Standard Combiner:
cmb_fast = imc.Combiner(stack)
out_sc_fast = cmb_fast.combine(
"median",
rejectors=imc.SigClip(sigma=(3.0, 3.0), maxiters=5, cenfunc="median"),
diagnostics=None,
)
print(f"Production output shape: {out_sc_fast.shape}")Production output shape: (8, 8)The short Chained Combiner snippets below exist only to inspect masks, clipping bounds, iteration counts, and status flags.
How rejection is evaluated
Every rejection method is evaluated independently at each output position. For one position, call the values from the input images a column. Input-masked and non-finite values are excluded before the algorithm begins. The equations below describe one such column; the same calculation runs for every output pixel.
2. Sigma Clipping
Iteratively clips values whose residual from the per-pixel center exceeds sigma * spread. Here sigma is the user-supplied multiplier, such as 3.0 (not the measured standard deviation). The measured spread comes from stdfunc (e.g., standard deviation or MAD).
For the detailed explanation (iterative equations, clip_cen/stdfunc roles, and rollback rules), see the SigmaClip algorithm.
Best for: general-purpose rejection with near-Gaussian noise.
c_sc = imc.Combiner(stack).reject(
imc.SigClip(sigma=(3.0, 3.0), maxiters=5, cenfunc="median")
)
out_sc = c_sc.combine("median")
mask_sc, std_sc, low_sc, upp_sc, nit_sc, code_sc = c_sc.last_reject
print(f"SigClip - rejected {mask_sc.sum():4d} pixels, max iters used: {nit_sc.max()}")SigClip - rejected 70 pixels, max iters used: 3low_sc and upp_sc are the lowest and highest retained finite sample values after clipping. The clipping thresholds instead include the center: center - sigma[0] * spread and center + sigma[1] * spread. Use low_sc and upp_sc to inspect the range of surviving values.
The nit array shows how many clipping passes were needed per output element:
unique_iters, counts = np.unique(nit_sc, return_counts=True)
for k, v in zip(unique_iters, counts):
print(f" nit={k}: {v} pixels") nit=2: 58 pixels
nit=3: 6 pixels3. CCD Clipping
Uses a CCD noise calculation instead of measuring the spread from the stack. The threshold changes with the local count level. This example has no zero or scale adjustment, so stack is in input DN. See the CCDClip algorithm for the detailed explanation.
Best for: data with well characterized CCD parameters.
c_cc = imc.Combiner(stack).reject(
imc.CcdClip(
sigma=(3.0, 3.0),
gain=1.0,
rdnoise=5.0,
maxiters=5,
)
)
out_cc = c_cc.combine("median")
mask_cc, std_cc, low_cc, upp_cc, nit_cc, code_cc = c_cc.last_reject
print(f"CcdClip - rejected {mask_cc.sum():4d} pixels, max iters used: {nit_cc.max()}")CcdClip - rejected 62 pixels, max iters used: 2low_cc and upp_cc are the lowest and highest retained finite sample values after clipping. They are sample extrema, not the thresholds used to reject values. CCD thresholds depend on the clipping center and each image’s noise; multiplying the reported std_cc by sigma does not reconstruct them when per-image noise differs.
zero= and scale= normalize the data before every rejection method. They are not CCDClip options. The direct CCDClip kernel has separate scales and zeros arguments because it needs them to calculate each image’s CCD noise. See Zero + Scale: CCDClip and direct kernels for their meaning and the sigscale rule.
The std diagnostic reports the unscaled reference noise evaluated at the clip_cen center. The rejection residuals use the cenfunc center and per-image noise. A non-finite reference variance stops iteration; ddof does not affect the CCD noise formula. This diagnostic is not the uncertainty of the combined image.
4. MinMax Clipping
Rejects the smallest and largest unmasked values at each output element. When n_min or n_max is greater than or equal to 1, it is interpreted as a frame count. When it is between 0 and 1, it is interpreted as a fraction of the original stack size and converted with IRAF-style truncation, int(N * fraction + 0.001). It makes a single rank-selection pass. The effective count is reduced proportionally when some input samples are masked or non-finite; see the equation below.
Best for: fixed N where you know exactly how many frames to trim, or as a first pass to knock out the worst extremes before statistical clipping.
See the MinMax algorithm for the exact per-pixel rank count calculation when a column has masked or non-finite values.
c_mm = imc.Combiner(stack).reject(imc.MinMaxClip(n_min=1, n_max=1))
out_mm = c_mm.combine("median")
mask_mm, std_mm, low_mm, upp_mm, nit_mm, code_mm = c_mm.last_reject
print(f"MinMaxClip - rejected {mask_mm.sum():4d} pixels")MinMaxClip - rejected 128 pixelslow_mm and upp_mm are inclusive retained-value bounds. Values equal to a bound are kept.
For a 25-frame stack, n_max=0.08 rejects int(25 * 0.08 + 0.001) == 2 maximum-side samples per output element:
c_mm_frac = imc.Combiner(stack).reject(imc.MinMaxClip(n_min=0, n_max=0.08))
print(f"Fractional MinMaxClip rejected {c_mm_frac.mask_rej.sum()} pixels")Fractional MinMaxClip rejected 128 pixels5. Linear Clipping
LinearClip applies a configurable affine interval around the current center. The default LinearClip() explicitly disables clipping. The 1-D kernel has a different default behavior, described below. maxiters=1 makes a single pass. See the LinearClip algorithm for the retained interval and iteration rules.
c_lin = imc.Combiner(stack).reject(
imc.LinearClip(
low_scale=1.0,
low=25.0,
upp_scale=1.0,
upp=25.0,
maxiters=2,
)
)
out_lin = c_lin.combine("median")
mask_lin, std_lin, low_lin, upp_lin, nit_lin, code_lin = c_lin.last_reject
print(f"LinearClip - rejected {mask_lin.sum():4d} pixels")LinearClip - rejected 67 pixelslow_lin and upp_lin are the final computed retained bounds. Unlike SigClip and CcdClip, this method does not estimate a standard deviation; it only applies the configured affine envelope around the chosen center. Use cenfunc="mean", "median", or "lmedian"/"lmed", and use nkeep, maxrej, and revert_on_nkeep when the same rollback safeguards are needed.
6. Percentile Clipping
IRAF-style pclip is a single-pass rejection method. It sorts each stack per-output stack, converts the scalar frac value to a rank offset around the median, uses that ranked sample to estimate a spread, then applies the lower and upper sigma thresholds.
Best for: IRAF compatibility or fast, single-pass clipping when a sorted rank offset is a good enough spread proxy for the stack.
See the PClip algorithm for rank resolution, the spread estimate, endpoint thresholds, and nkeep restoration.
c_pc = imc.Combiner(stack).reject(
imc.PClip(frac=-0.5, sigma=(3.0, 3.0))
)
out_pc = c_pc.combine("median")
mask_pc, std_pc, low_pc, upp_pc, nit_pc, code_pc = c_pc.last_reject
print(f"PClip - rejected {mask_pc.sum():4d} pixels")PClip - rejected 246 pixelslow_pc and upp_pc are the lowest and highest retained finite sample values after pclip rejection. Use mask_pc to identify which samples survived.
7. Algorithms
Every method below operates independently on one output column: the values at one spatial position across input images. Input-masked and non-finite values are excluded before the method begins.
Where the methods differ from common textbook versions, they follow IRAF imcombine; see its documentation.
SigmaClip
At iteration \(t\), let \(n_t\) be the number of retained finite, unmasked values, \(c_t\) the cenfunc (for example, mean or median) center used for clipping, and \(c_{\mathrm{clipcen},t}\) the center used only to estimate the spread. clip_cen selects the function for \(c_{\mathrm{clipcen},t}\); when clip_cen=None (the default), \(c_{\mathrm{clipcen},t}=c_t\). stdfunc selects the spread \(d_t\) (standard deviation or MAD).
\[ c_t - {\rm lsigma} \times d_t \mathrel{\leq} x \mathrel{\leq} c_t + {\rm hsigma} \times d_t. \]
Here \({\rm lsigma}\) and \({\rm hsigma}\) are the lower and upper sigma multipliers. Values equal to either threshold are retained. All values outside the interval are tentatively rejected together, then the center and spread are recomputed from the survivors.
For stdfunc="std", the spread is the standard deviation about clip_cen:
\[ d_t = \sqrt{\frac{\sum_{\text{retained }x}(x-c_{\mathrm{clipcen},t})^2} {n_t-\mathrm{ddof}}}. \]
For stdfunc="mad", it is
\[ d_t = 1.4826\,\operatorname{median}_{\text{retained }x} \left|x-c_{\mathrm{clipcen},t}\right|. \]
ddof applies only to the standard-deviation form; it is ignored for MAD. Iteration stops when no new values are rejected or maxiters is reached.
CCDClip
CCDClip calculates noise in input DN. Let \(C\) be the current cenfunc center of the normalized values. For image \(i\), it first converts \(C\) to an input-DN signal \(S_i\). It then converts that image’s noise standard deviation back to normalized units:
\[ S_i = \max(0, s_i(C + z_i)), \qquad \sigma_i = \frac{\sqrt{\max((r/g)^2, 10^4/\mathrm{MAX\_REAL}) + S_i/g + (qS_i)^2}}{s_i}. \]
Here MAX_REAL = 3.402823466e38 is IRAF’s single-precision limit. s_i and z_i are values used only by the direct CCDClip kernel. Do not pass public scale= and zero= values as scales and zeros. The exact mapping is in Zero + Scale: CCDClip and direct kernels. g is gain, r is rdnoise, and q is snoise: the fractional sensitivity noise. For example, flat-field uncertainty adds the \((qS_i)^2\) term to the variance. If no per-image vectors are active, CCDClip uses \(S=\max(0,C)\).
At each iteration, it then does the following:
Sort retained pairs \((x_i, i)\) by value, keeping each sample’s original image index.
Scan from the lowest value inward. Reject while
\[R_{\mathrm{low},i} = \frac{C-x_i}{\sigma_i} > \lambda_-,\]
and stop this side at the first retained sample.
Independently scan from the highest value inward. Reject while
\[R_{\mathrm{high},i} = \frac{x_i-C}{\sigma_i} > \lambda_+,\]
and stop this side at the first retained sample.
Apply the shared iterative safeguards. If the pass is accepted, recompute \(C\) from the survivors.
Without per-image vectors, every sample has the same CCD sigma. With them, each sample has its own \(\sigma_i\). This is why CCDClip cannot test every sample independently. For example, suppose \(C=100\) and the lower multiplier is 3:
| sorted value | sample sigma | lower residual |
|---|---|---|
| 50 | 20 | 2.5 (retained) |
| 60 | 10 | 4.0 (would fail if tested) |
IRAF sees 50 first, retains it, and stops the low-side scan. It never tests 60, so 60 is retained too. The noisier, more extreme value shields the nearer value. Testing every sample independently would reject 60 and would not reproduce IRAF CCDClip.
MinMaxClip
MinMax is a single sorted-rank pass. Let \(N\) be the original number of input images and \(M\) the number of finite, unmasked values at this output pixel. For a requested n_min or n_max value \(a\), the public API first resolves a frame count
\[ A = \begin{cases} \operatorname{int}(N a + 0.001), & 0 \leq a < 1,\\ \operatorname{int}(a), & a \geq 1. \end{cases} \]
The kernel then rejects \(\operatorname{int}(M A/N + 0.001)\) values from that tail. Thus an input mask or non-finite value can reduce the number rejected at a particular output pixel. It sorts the usable values once, rejects the requested low-end ranks, then rejects the requested high-end ranks from the remaining values. There is no center, spread estimate, or iteration.
LinearClip
At iteration \(t\), compute the selected center \(c_t\) from the current survivors and form
\[ L_t = a_-c_t-b_-, \qquad U_t = a_+c_t+b_+, \]
where \(a_-\), \(b_-\), \(a_+\), and \(b_+\) are low_scale, low, upp_scale, and upp. Retain exactly the values satisfying \(L_t \leq x \leq U_t\); equality is retained. All out-of-range values are tentatively rejected together. An accepted pass recomputes the center from the survivors; a pass violating the shared safeguards is undone in full. The default maxiters=1 performs one clipping pass.
The stack wrapper linearclip() and LinearClip() explicitly bypass this equation when (low_scale, low, upp_scale, upp) = (1, 0, 1, 0): they return an all-false rejection mask, NaN bounds, and zero iterations. This disables clipping; it is not the mathematical result of substituting those parameters. linearclip_1d() instead applies the equation, so the default interval is [c, c]. It keeps only samples equal to c, unless nkeep or maxrej causes the pass to be reverted. Specify non-default bounds explicitly when comparing the stack and 1-D APIs.
PClip
PClip is a sorted-rank method that estimates a spread from one ranked sample, rather than from all deviations. Let \(N\) be the original image count. First resolve the signed rank offset \(k\): if \(|\mathrm{frac}|<1\), multiply it by \(\max(\lfloor N/2\rfloor,1)\), truncate toward zero, then clamp it to the nonzero interval from \(-\max(\lfloor N/2\rfloor,1)\) through \(+\max(\lfloor N/2\rfloor,1)\). Larger absolute frac values are truncated and clamped directly.
For the \(M\) finite, unmasked values, sort \(x_{(1)} \leq \cdots \leq x_{(M)}\) and compute their usual median \(m\) (the average of the two middle values for even \(M\)). In one-based indexing, let \(b=1+\lfloor M/2\rfloor\). The selected rank is
\[ r = \operatorname{clip}\!\left( \begin{cases} b-1+k, & M\text{ even and }k<0,\\ b+k, & \text{otherwise}, \end{cases}, 1, M\right). \]
Its spread estimate is
\[ d = \operatorname{sign}(k)\,[x_{(r)}-m]. \]
It rejects low and high sorted endpoints while, respectively,
\[ \frac{m-x_{(j)}}{d} \geq \lambda_- \quad\text{or}\quad \frac{x_{(j)}-m}{d} \geq \lambda_+. \]
Unlike SigmaClip and CCDClip, equality is rejected here. The method does not iterate. If this would leave fewer than nkeep samples, it restores the least-extreme rejected endpoints using IRAF’s rank-based restoration rule. With fewer than three usable samples, fewer than nkeep + 1 usable samples, or a zero rank-based spread, PClip makes no rejection.
Grow rejected samples (all methods)
grow expands the rejection mask after any rejection decision: SigClip, CcdClip, LinearClip, MinMaxClip, or PClip. Pass it either on the rejector or to Combiner.reject(..., grow=...). It is useful when a cosmic ray or bad pixel should also remove nearby samples in the same input frame. The expansion is spatial only: axis 0 is the stack axis and is not grown across. The 1-D kernels do not accept grow, because they have no spatial neighbors.
The radius is exact Euclidean distance on the integer pixel grid. In 2-D, grow=1 reaches the four direct neighbors, while grow=np.sqrt(2) also reaches the diagonal neighbors. grow=None is the default and skips this extra calculation.
toy = np.ones((10, 5, 5), dtype="float32")
toy[0, 2, 2] = 100.0
c_grow = imc.Combiner(toy).reject(
imc.SigClip(sigma=3.0, maxiters=5, ddof=0, nkeep=0),
grow=1,
)
print(c_grow.mask_rej[0].astype(int))
print(f"grow-added output flags: {np.count_nonzero(c_grow.output_flags & imc.OutputFlags.GROW)}")[[0 0 0 0 0]
[0 0 1 0 0]
[0 1 1 1 0]
[0 0 1 0 0]
[0 0 0 0 0]]
grow-added output flags: 4mask_rej includes the algorithm’s exclusions and samples added by growth. Sigma, CCD, MinMax, and PClip masks also include prior masks and non-finite samples, so those entries can seed growth too. LinearClip excludes input-mask entries from its returned mask; its disabled default returns an empty mask. Later combine calls exclude every value in the running mask. low, upp, nit, and std describe the rejection before growth; std is None for methods without a spread model. imc.OutputFlags.GROW (value 16) marks an output element where growth added at least one rejected sample.
The cost of grow depends mostly on how many samples were rejected and how large the growth footprint is. The following local timing example uses SigmaClip merely as one representative workload; growth has the same post-rejection role for every method. Treat the values as machine-local guidance, not a fixed benchmark.
bench_rng = np.random.default_rng(20250311)
bench = bench_rng.normal(1000.0, 20.0, (31, 256, 256)).astype("float32")
bench[0, ::32, ::32] += 800.0
bench[1, 16::32, 16::32] -= 600.0
timings = [
bench_case("median only", median_only),
bench_case("sigclip median, grow=None", lambda: sigclip_median(None)),
bench_case("sigclip median, grow=1", lambda: sigclip_median(1)),
bench_case("sigclip median, grow=sqrt(2)", lambda: sigclip_median(np.sqrt(2))),
bench_case("sigclip median, grow=2", lambda: sigclip_median(2)),
]
baseline_ms = timings[1][1]
print("| case | rejected samples | median time | relative to no grow |")
print("|---|---:|---:|---:|")
for label, ms, (_, mask) in timings:
n_rej = "-" if mask is None else str(int(mask.sum()))
rel = "-" if label == "median only" else f"{ms / baseline_ms:.2f}x"
print(f"| {label} | {n_rej} | {ms:.2f} ms | {rel} |")
assert timings[0][2][0].shape == bench.shape[1:]
assert np.isfinite(timings[1][2][0]).all()
assert timings[2][2][1].sum() >= timings[1][2][1].sum()| case | rejected samples | median time | relative to no grow |
|---|---:|---:|---:|
| median only | - | 7.85 ms | - |
| sigclip median, grow=None | 3898 | 36.28 ms | 1.00x |
| sigclip median, grow=1 | 19350 | 45.36 ms | 1.25x |
| sigclip median, grow=sqrt(2) | 34618 | 44.95 ms | 1.24x |
| sigclip median, grow=2 | 49804 | 38.52 ms | 1.06x |Exact timings vary with the operation, dtype, image count, image shape, rejection density, grow radius, CPU, and thread settings. Larger radii and dense rejection masks cost more because more neighboring samples must be marked. Use grow when neighboring samples should be removed with the rejected sample; too large a radius can discard useful data around real sources or dense artifacts.
8. Compare Outputs
Each Chained Combiner example stores rejected values in its mask. Its later combine call therefore uses only finite, unmasked, unrejected values:
cmb_raw = imc.Combiner(stack)
out_raw = cmb_raw.combine("median")
print(f"True value: {true_val:.1f}")
print(f"No rejection: mean={out_raw.mean():.2f}, std={out_raw.std():.2f}")
print(f"SigClip (3): mean={out_sc.mean():.2f}, std={out_sc.std():.2f}")
print(f"CcdClip (3): mean={out_cc.mean():.2f}, std={out_cc.std():.2f}")
print(f"LinearClip: mean={out_lin.mean():.2f}, std={out_lin.std():.2f}")
print(f"MinMaxClip (1,1): mean={out_mm.mean():.2f}, std={out_mm.std():.2f}")
print(f"PClip (-0.5): mean={out_pc.mean():.2f}, std={out_pc.std():.2f}")True value: 100.0
No rejection: mean=100.06, std=1.32
SigClip (3): mean=99.80, std=1.30
CcdClip (3): mean=99.82, std=1.30
LinearClip: mean=99.81, std=1.30
MinMaxClip (1,1): mean=100.06, std=1.32
PClip (-0.5): mean=99.81, std=1.449. Chain Rejectors
All rejection algorithms are frozen-dataclass Rejector subclasses. For production reductions, pass a list of rejectors to Standard Combiner:
# Trim the extremes first, then sigma-clip the rest
cmb_combo = imc.Combiner(stack)
out_combo = cmb_combo.combine(
"median",
rejectors=[
imc.MinMaxClip(n_min=1, n_max=1),
imc.SigClip(sigma=3.0, maxiters=5),
],
diagnostics=None,
)
print(f"MinMaxClip + SigClip: mean={out_combo.mean():.2f}")MinMaxClip + SigClip: mean=100.03Use Chained Combiner for the same sequence when you need intermediate masks and diagnostics retained on the combiner.
10. Variance Maps
Variance is computed from the final valid values, not from the rejection threshold model. In NumPy terms:
arr_eff = np.where(mask | rejected, np.nan, stack)
var = np.nanvar(arr_eff, axis=0, ddof=ddof)With Combiner, the current mask already includes input masks and all rejections:
var, mean = c_sc.variance(ddof=1, return_mean=True)
err = np.sqrt(var)
print(
f"mean={np.nanmean(mean):.2f}, "
f"variance mean={np.nanmean(var):.2f}, "
f"error mean={np.nanmean(err):.2f}"
)mean=99.91, variance mean=24.41, error mean=4.89Use return_mean=True when you need both variance and mean from the same final valid values. np.sqrt(var) measures their scatter; it is not the uncertainty of the combined image.
11. lmedian
For an even number of samples, median averages the middle two values. lmedian takes the lower middle value, following IRAF. You can also use lmedian or lmed for cenfunc and clip_cen; that changes rejection thresholds. combine="lmedian" instead changes the final stack reduction.
Because lmedian selects an input sample, it can keep an integer dtype. This applies to ndcombine(..., combine="lmedian") or "lmed" only when there is no mask, rejection, zero, scale, or weight. For an even-sized uint16 stack, ordinary median may produce a value ending in .5 and therefore needs a floating output. Lower median can stay uint16 and use less storage. Calls with masks, rejection, zero/scale, or Combiner use the normal floating workspace.
even_stack = stack[:10] # N=10 (even)
cmb_even = imc.Combiner(even_stack)
out_med = cmb_even.combine("median")
out_lmed = imc.ndcombine(even_stack, combine="lmedian")
diff = out_lmed - out_med
print(f"lmedian - median: min={diff.min():.4f}, max={diff.max():.4f}")
print("(lmedian <= median everywhere for even N)")lmedian - median: min=-2.2517, max=-0.0008
(lmedian <= median everywhere for even N)The combine string names and aliases are the same in Combiner and ndcombine:
| Method | Also accepts |
|---|---|
mean |
average, avg |
median |
med |
lmedian |
lmed |
sum |
|
min |
|
max |
|
variance |
var |
| weighted mean | use combine="mean", weight=w |
12. Performance Tips: direct kernel calls
Direct kernel calls are intended for already-prepared arrays and custom high-throughput layers. They do not keep Combiner state, so you must apply rejection masks yourself:
mask_direct, std_direct, low_direct, upp_direct, nit_direct, code_direct = imck.sigclip(
stack, sigma=(3.0, 3.0), maxiters=5, cenfunc="median"
)
out_direct = imck.sigclip_combine(
stack,
combine="median",
sigma=(3.0, 3.0),
maxiters=5,
cenfunc="median",
)
np.testing.assert_allclose(out_direct, out_sc)
print("direct kernel calls match the Chained Combiner result.")direct kernel calls match the Chained Combiner result.Standard Combiner should remain the default API unless profiling shows that wrapper overhead matters.
13. Quick Reference
| Algorithm | Deterministic? | Adapts to noise? | Needs parameters |
|---|---|---|---|
SigClip |
No | Yes (empirical sigma) | sigma, maxiters, cenfunc |
CcdClip |
No | Yes (CCD model) | sigma, gain, rdnoise |
LinearClip |
Yes | Recomputed center, no spread estimate | low_scale, low, upp_scale, upp, maxiters, cenfunc |
MinMaxClip |
Yes | No | n_min, n_max |
PClip |
Yes | Approximate spread from sorted rank | frac, sigma, nkeep |