04 Zero Scale

zero= subtracts one offset from each input frame. scale= then divides each frame by one factor. Together they calculate

\[ \mathrm{normalized} = \frac{\mathrm{inputdata} - \mathrm{zero\_offset}} {\mathrm{scale\_factor}}. \]

Standard Combiner and ndcombine do this before rejectors passed to the same combine call. Chained Combiner does it when you call zero_scale(...). This tutorial uses Standard Combiner examples. The zero and scale arguments accept three forms:

Form	Example	When to use
Numeric array	scale=illumination	You have the values from a separate measurement
String statistic	scale="median", "med_sc"	Auto-compute a per-frame statistic over all pixels
Callable	scale=lambda p: np.percentile(p, 10)	Custom estimator (e.g. robust background)

By default, each zero and scale value is expressed relative to the zeroth frame (zero_to_0th=True, scale_to_0th=True). If the values you provide are \(z_i\) and \(s_i\), the values actually used are

\[ \mathrm{zero\_offset}_i = z_i-z_0, \qquad \mathrm{scale\_factor}_i = s_i/s_0. \]

The zeroth frame is therefore unchanged. Set either option to False when you want to use the values exactly as provided.

This is imcombiners’ normalization convention. IRAF’s joint zero/scale calculation can differ when both are derived from statistics and scales vary; do not assume full normalization parity from the reference-frame options.

imc.resolve_zero_scale(...) can show the one zero offset and scale factor that will be used for each frame. Usually, pass zero= and scale= directly to Standard Combiner instead.

import numpy as np
import imcombiners as imc

1. Synthetic Setup

This example simulates flat-field exposures taken under different lamp brightnesses. Each frame has 2% pixel-to-pixel variation, and its overall level differs from the other frames.

rng = np.random.default_rng(20250311)
n_frames, ny, nx = 10, 32, 32

# Spatially non-uniform sensitivity pattern (same for all frames)
spatial_pattern = rng.normal(1.0, 0.02, (ny, nx)).astype("float32")

# Each frame has a different overall illumination level
illumination = np.array(
    [1000, 1200, 950, 1100, 1050, 1300, 1150, 900, 1080, 1120],
    dtype="float32",
)

stack = np.stack([
    spatial_pattern * ill + rng.normal(0, 5.0, (ny, nx)).astype("float32")
    for ill in illumination
])

print("Per-frame means (dominated by illumination):")
print(stack.mean(axis=(1, 2)).round(0))

Per-frame means (dominated by illumination):
[ 999. 1200.  950. 1099. 1050. 1300. 1150.  899. 1079. 1120.]

Without normalization, frames with different illumination would be combined unequally and the result would be biased toward the brightest frames:

cmb_raw = imc.Combiner(stack)
out_raw = cmb_raw.combine("median")
print(f"Raw median: mean={out_raw.mean():.1f}  (expect ~illumination.mean()={illumination.mean():.0f})")

Raw median: mean=1089.4  (expect ~illumination.mean()=1085)

The goal is to normalize each frame to unit illumination before combining.

2. Numeric Values

Use this when you have the per-frame levels from an independent measurement (e.g. flat frames with different illumination levels):

cmb_num = imc.Combiner(stack)
out_num = cmb_num.combine(
    "median",
    scale=illumination,
    scale_to_0th=False,
)
print(f"Numeric scale: mean={out_num.mean():.4f}  (expect ~1.0)")

Numeric scale: mean=0.9996  (expect ~1.0)

The examples in this section use scale_to_0th=False to make the combined image have unit normalization. With the default scale_to_0th=True, the first scale value is the reference, so the output stays on the zeroth frame’s scale. This may be of interest if you want to preserve the Poisson noise level to a reasonable degree. Similarly, use zero_to_0th=False when zero is an absolute sky background that should be subtracted in full.

3. String Statistics

Pass one of "mean", "median", "sum", "min", or "max" as zero= or scale=, and Standard Combiner computes that statistic over each plane’s pixels automatically. Sigma-clipped statistics are available as "sigclip_mean" / "mean_sc" and "sigclip_median" / "median_sc" / "med_sc".

Built-in statistics use finite samples only: both NaN and infinity are excluded. For example, a plane containing [1, inf, 3] gives a mean scale of 2 with scale_to_0th=False. Custom callables choose their own nonfinite policy.

The sigma-clipped aliases use sigma=3, maxiters=5, cenfunc="median", clip_cen=None, and ddof=0 by default. Tune them with zero_sigclip_kwargs or scale_sigclip_kwargs; cenfunc and clip_cen also accept lower median as "lmedian"/"lmed".

These aliases clip float64 samples, with no minimum-survivor rollback. Rejections are cumulative: a rejected value stays excluded if later thresholds move. Zero spread follows the squared threshold rule instead of stopping by itself; floating-point underflow still limits comparisons. If no samples survive, the statistic is NaN and normal zero/scale validation rejects it.

# Divide each frame by its per-plane median — robust to pixel outliers
cmb_str = imc.Combiner(stack)
out_str = cmb_str.combine(
    "median",
    scale="median",
    scale_to_0th=False,
)
print(f"scale='median': mean={out_str.mean():.4f}  (expect ~1.0)")

scale='median': mean=1.0001  (expect ~1.0)

cmb_sc = imc.Combiner(stack)
out_sc = cmb_sc.combine(
    "median",
    scale="med_sc",
    scale_to_0th=False,
    scale_sigclip_kwargs={"sigma": (2.5, 4.0), "maxiters": 3},
)
print(f"scale='med_sc': mean={out_sc.mean():.4f}")

scale='med_sc': mean=0.9999

For this finite synthetic stack, this is equivalent to computing the medians manually and passing them:

per_plane_median = np.median(stack, axis=(1, 2))
cmb_equiv = imc.Combiner(stack)
out_equiv = cmb_equiv.combine(
    "median",
    scale=per_plane_median,
    scale_to_0th=False,
)
np.testing.assert_allclose(out_str, out_equiv, rtol=1e-5)
print("Equivalent to manual per-plane median.")

Equivalent to manual per-plane median.

4. Callable Statistics

For a custom estimator — anything that maps one input plane with shape (*spatial) to a scalar:

# Robust level estimate: use the central 80th-percentile range midpoint
def robust_level(plane):
    lo, hi = np.percentile(plane, [10, 90])
    return float(0.5 * (lo + hi))

cmb_callable = imc.Combiner(stack)
out_callable = cmb_callable.combine(
    "median",
    scale=robust_level,
    scale_to_0th=False,
)
print(f"Callable scale: mean={out_callable.mean():.4f}  (expect ~1.0)")

Callable scale: mean=1.0002  (expect ~1.0)

5. Background Subtraction

zero subtracts a per-frame offset before dividing by scale. A typical use is removing per-frame sky backgrounds from science frames.

# Simulate science frames: source + sky background varying per frame
sky = np.array([200, 210, 195, 205, 200, 215, 208, 192, 203, 207], dtype="float32")
source = 500.0

sci_stack = np.stack([
    rng.normal(source, 15.0, (ny, nx)).astype("float32") + s
    for s in sky
])

# Subtract the known per-frame sky, then combine
cmb_zero = imc.Combiner(sci_stack)
out_zero = cmb_zero.combine(
    "mean",
    zero=sky,
    zero_to_0th=False,
)
print(f"After sky subtraction: mean={out_zero.mean():.1f}  (expect ~{source:.0f})")

After sky subtraction: mean=500.1  (expect ~500)

6. Zero + Scale

Apply sky subtraction and flux normalization together in one Standard Combiner call. The order is always \((\mathrm{inputdata}-\mathrm{zero})/\mathrm{scale}\):

cmb_both = imc.Combiner(stack)
out_both = cmb_both.combine(
    "median",
    zero="median",
    scale="mean",
    rejectors=imc.SigClip(sigma=3.0, maxiters=5),
    diagnostics=None,
)
print(f"zero + scale + sigclip median: mean={out_both.mean():.4f}")

zero + scale + sigclip median: mean=916.6476

7. Practical Notes

• In Standard Combiner, zero and scale affect that call only. With diagnostics=None, they do not change cmb.arr.
• Every scale value must be finite and nonzero. Otherwise, the call raises ValueError.
• Use Chained Combiner (cmb.zero_scale(...).reject(...).combine(...)) to inspect normalized data before combining.

8. CCDClip and direct normalized kernels

CCDClip clips normalized values, but it calculates CCD noise in input DN. When the sigscale rule enables the per-image correction, CCDClip converts the current normalized center back to input DN separately for each image before it calculates that image’s sigma. It uses this conversion only for the sigma calculation. It does not normalize the data again.

Let zero_offset_i and scale_factor_i be the public values after any zero_to_0th and scale_to_0th rebasing. The public normalized data are

\[ \mathrm{normalized}_i = \frac{\mathrm{inputdata}_i - \mathrm{zero\_offset}_i} {\mathrm{scale\_factor}_i}. \]

The direct kernels.ccdclip API uses scales and zeros for that sigma calculation. They are not the same values as public scale= and zero=:

\[ s_i = \mathrm{scale\_factor}_i, \qquad z_i = \frac{\mathrm{zero\_offset}_i}{\mathrm{scale\_factor}_i}, \qquad \mathrm{inputdata}_i = s_i(\mathrm{normalized}_i + z_i). \]

Usually, pass public zero= and scale= to Combiner.combine(...) or ndcombine(...); do not supply these direct-kernel vectors. If you call a CCDClip kernel on data that you already normalized, supply them explicitly:

normalized = np.ones((3, 1, 1), dtype="float32")
mask_ccd, *_ = imc.kernels.ccdclip(
    normalized,
    scales=np.array([1.0, 1.5, 0.8]),
    zeros=np.array([0.0, 4.0, -2.0]),
    gain=2.0,
    rdnoise=5.0,
)

IRAF uses the per-image sigma calculation only when both conditions are true:

• The resolved scale factors are not all equal.
• At least one scale factor differs from \(1\) by more than sigscale (default 0.1).

ndcombine(..., sigscale=...) exposes this rule. Standard Combiner uses the IRAF default. When either condition is false, CCDClip still clips the normalized data, but calculates sigma as if no per-image scale correction were active.

9. Prepared Workspace Check

For simple fixed normalization, a prepared workspace can reproduce the same result:

manual = np.ascontiguousarray(stack / illumination[:, None, None], dtype=np.float32)
out_low = imc.Combiner(manual, copy=False, validate=False).combine("median")

np.testing.assert_allclose(out_num, out_low, rtol=1e-5)
print("prepared workspace matches the numeric-scale result.")

prepared workspace matches the numeric-scale result.