07 IRAF Compatibility

This tutorial translates IRAF-style combine settings into Standard Combiner. Use compact ndcombine() wrapper when a compact IRAF-like function call is useful, Chained Combiner when diagnostics are part of the task, and direct kernel calls only for custom performance-sensitive layers.

import numpy as np
import imcombiners as imc
import imcombiners.kernels as imck

1. Standard Combiner

rng = np.random.default_rng(20250311)
stack = rng.normal(100, 10, (15, 32, 32)).astype("float32")

cmb = imc.Combiner(stack)
out = cmb.combine(
    "median",
    rejectors=imc.SigClip(
        sigma=(3.0, 3.0),
        maxiters=5,
        cenfunc="median",
        nkeep=1,
        revert_on_nkeep=True,
    ),
    diagnostics=None,
)
print(out.shape, out.dtype)

(32, 32) float32

The equivalent compact ndcombine() wrapper spelling keeps the IRAF-style shape:

out_one = imc.ndcombine(
    stack,
    combine="median",
    reject="sigclip",
    sigma=(3.0, 3.0),
    maxiters=5,
    cenfunc="median",
    nkeep=1,
    revert_on_nkeep=True,
)
np.testing.assert_allclose(out, out_one)
print("Standard Combiner and compact ndcombine() wrapper agree.")

Standard Combiner and compact ndcombine() wrapper agree.

reject accepts the same strings: "sigclip", "ccdclip", "minmax", "pclip".

For dtype-sensitive integer stacks, pure combine="lmedian" or combine="lmed" preserves the input dtype when no mask, rejection, zero, scale, or weight is used. Once any of those workflow features are enabled, ndcombine uses the normal floating workspace.

Accepted public input dtypes are uint8, uint16, int16, int32, float32, and float64. Other dtypes are not silently downgraded. If your data are int64, uint64, float128, or another unsupported dtype, cast them explicitly before calling ndcombine. There is no full-precision int64 combine path because the package does not provide a float128 workspace.

2. Kwarg Map

IRAF `imcombine` parameter/value	`ndcombine(...)` bridge kwarg	`Combiner` equivalent
`lthreshold=lo`, `hthreshold=hi`	`thresholds=(lo, hi)`	`.threshold(lo, hi)`
`zero=z`, `scale=s`	`zero=z, scale=s`	`.zero_scale(zero=z, scale=s)`
`combine="average"`	`combine="mean"` or `"average"`	`.combine("mean")`
`combine="median"`	`combine="median"`	`.combine("median")`
`combine="sum"`	`combine="sum"`	`.combine("sum")`
`reject="none"`	`reject=None`	omit `.reject(...)`
`reject="sigclip"`, `lsigma=lo`, `hsigma=hi`	`reject="sigclip", sigma=(lo, hi)`	`.reject(SigClip(sigma=(lo, hi)))`
`reject="avsigclip"`	no direct equivalent	`SigClip(cenfunc="mean")` does not implement IRAF’s line-based noise estimate
`reject="ccdclip"`, `rdnoise=r`, `gain=g`, `snoise=s`	`reject="ccdclip", rdnoise=r, gain=g, snoise=s`	`.reject(CcdClip(rdnoise=r, gain=g, snoise=s))`
`nlow=lo`, `nhigh=hi` with `reject="minmax"`	`reject="minmax", n_minmax=(lo, hi)`	`.reject(MinMaxClip(n_min=lo, n_max=hi))`
`pclip=p` with `reject="pclip"`	`reject="pclip", pclip=p`	`.reject(PClip(frac=p))`
`nkeep=n`	`nkeep=n`	`SigClip(..., nkeep=n)` or `CcdClip(..., nkeep=n)`
`mclip=yes` / `mclip=no`	`cenfunc="median"` / `cenfunc="mean"`	`SigClip(cenfunc=...)` or `CcdClip(cenfunc=...)`; `cenfunc="lmedian"`/`"lmed"` is also available when lower-median centering is wanted
`niterate=n`	`maxiters=n`	`SigClip(maxiters=n)` or `CcdClip(maxiters=n)`
`grow=r`	`grow=r`	`.reject(SigClip(...), grow=r)` or `SigClip(..., grow=r)`
`weight=w`, `combine="average"`	`weight=w, combine="mean"`	`.combine("mean", weight=w)`
`rejmasks`, `nrejmasks`, `expmasks`, `sigmas` diagnostic products	`diagnostics="simple"`	returns `(combined, mask_rej, mask_thresh, std, low, upp, nit, output_flags)`, not IRAF-format files
`blank=value`	no `blank` keyword; fill empty output pixels in the caller	output-fill policy, not an input mask
`offsets=...`	prepare an offset-padded stack first	see the offset stacking tutorial, then combine the padded stack
`masktype`, `maskvalue`	prepare a boolean input mask first	pass that mask to `ndcombine(..., mask=mask)` or `Combiner(..., mask=mask)`
`sigscale`	`sigscale=0.1` for CCDClip	CCDClip applies IRAF’s scale-gated per-image noise vectors; scaled SigClip remains unsupported
`reject="crreject"`, `project`	no direct equivalent	handle upstream/downstream before `Combiner`

Thresholds are applied before zero/scale, rejection, and combination. Bounds are inclusive: values survive when lo <= value <= hi. With diagnostics="simple", mask_thresh records threshold exclusions. mask_rej also includes input masks, threshold exclusions, and non-finite samples for the supported rejection methods. Use diagnostics="full" and sample_flags to distinguish causes.

grow expands mask_rej spatially within each input frame after the rejection algorithm runs. The radius is exact Euclidean distance on the integer pixel grid: grow=1 reaches the four direct neighbors in 2-D, while grow=sqrt(2) also reaches diagonal neighbors. It does not grow across the stack axis. Because those methods include prior exclusions in mask_rej, input-mask, threshold, and non-finite entries can also seed growth; the separate input mask and mask_thresh arrays themselves are not modified. With diagnostics="simple", output_flags value 16 marks output elements where growth added at least one rejected sample.

IRAF’s blank sets the output where no input samples survive. It does not identify an input sentinel. If your combined result uses NaN for empty output pixels, replace those values in the caller. IRAF’s avsigclip estimates a noise proportionality factor across each output line; mean-centered SigClip is a different algorithm. See the IRAF imcombine documentation.

Convert IRAF mask files and masktype/maskvalue conventions to a boolean input mask before calling imcombiners. Prepare offsets upstream as shown in the offset stacking tutorial.

For CCDClip, ndcombine implements IRAF’s sigscale gate: it applies its resolved per-image zero/scale vectors to the CCD noise model only when scales differ from one by more than sigscale (default 0.1); sigscale=0 disables the correction. Scaled SigClip correction is unsupported.

The ndcombine(...) column maps IRAF options to a compact function call. The Combiner column shows Chained Combiner spelling because it maps each IRAF stage visibly; Standard Combiner can pass the same threshold and rejector specs directly to combine(...). IRAF compatibility is limited to the supported settings described here; bit-for-bit equivalence is not guaranteed for every mode. Sigma/CCD clipping recomputes the center during iterative rejection. Use revert_on_nkeep=True for IRAF-like nkeep rollback when matching IRAF-style rejection matters.

The default nkeep=1 is intentionally conservative: iterative clipping will not reject every finite sample in a per-output stack. Setting nkeep=0 removes that guard and can produce all-rejected output elements under aggressive clipping. With revert_on_nkeep=True, an iteration that would leave fewer than nkeep survivors is reverted in full for that pixel. Input-masked and non-finite samples remain excluded. Setting revert_on_nkeep=False makes clipping strict and may leave fewer than nkeep usable samples.

For median sigma/CCD clipping with no growth, IRAF computes the final median from the retained sorted window but skips moving a low-side rejected window before writing rejmasks. Therefore a low-side rejection can appear in IRAF’s rejection mask as the corresponding high sorted endpoint even though the final median used the low-side-clipped window. With diagnostics="simple", imcombiners marks the samples actually rejected by the clipping calculation, so diagnostic mask ordering can differ from IRAF.

For sigma clipping with clip_cen="mean", the spread uses a two-pass variance about the mean. Floating-point rounding near a clipping threshold can affect which samples survive.

For CCD images, common values are thresholds=(0, satlevel) to remove bad/cold pixels and saturated pixels, or thresholds=(0, 0.99 * satlevel) when near-saturated pixels should also be excluded.

3. IRAF Matching Settings

To reproduce supported IRAF-style output as closely as imcombiners can, translate the IRAF task parameters explicitly instead of relying on defaults:

out = imc.ndcombine(
    stack,
    combine="median",          # "mean"/"average", "median", "lmedian", or "sum"
    reject="sigclip",          # or "ccdclip", "minmax", "pclip", None
    sigma=(3.0, 3.0),          # IRAF lsigma, hsigma
    maxiters=5,                # imcombiners iteration limit
    cenfunc="median",          # IRAF mclip=yes; "mean" for mclip=no; "lmedian"/"lmed" for lower median
    clip_cen=None,             # None means use cenfunc each iteration
    nkeep=1,                   # IRAF nkeep
    revert_on_nkeep=True,      # revert an iteration that would violate nkeep
    zero_to_0th=True,          # reference-frame zero convention
    scale_to_0th=True,         # reference-frame scale convention
)

For combine="lmedian" with an even number of integer input images, pure ndcombine(stack, combine="lmedian") can preserve the integer dtype by taking the lower of the two middle samples. IRAF’s lmedian takes the lower sample only when exactly two values survive; it uses the usual median otherwise. Once masking, rejection, zero/scale, or weights are involved, imcombiners uses its normal floating workspace.

Use diagnostics=None when only the combined image is needed. Supported mean/median combinations with sigclip, ccdclip, minmax, or pclip avoid retaining a stack-sized rejection-mask result. Request diagnostics="simple" to inspect rejection products. Spatial growth (grow=...) requires a mask even when the final output is image-only.

pclip follows IRAF’s scalar parameter semantics. IRAF converts the task parameter into an integer rank offset from the sorted median, uses that sample’s distance from the median as a sigma estimate, then applies lsigma/hsigma clipping. In imcombiners, pass this value as pclip= in the compact ndcombine() wrapper, frac= in PClip, or frac= in imck.pclip. The value must be scalar.

IRAF-style sigscale correction is supported for CCDClip only; scaled SigClip correction is unsupported.

4. Performance Tips: direct kernel calls

Direct kernel calls expose the same Rust rejection primitives without the Combiner wrapper. They are useful when a pipeline has already prepared contiguous floating workspaces:

mask_low, *_ = imck.sigclip(
    stack,
    sigma=(3.0, 3.0),
    maxiters=5,
    cenfunc="median",
    nkeep=1,
    revert_on_nkeep=True,
)
out_low = imck.sigclip_combine(
    stack,
    combine="median",
    sigma=(3.0, 3.0),
    maxiters=5,
    cenfunc="median",
    nkeep=1,
    revert_on_nkeep=True,
)

np.testing.assert_allclose(out, out_low)
print("direct kernel calls reproduce the IRAF-style median.")

direct kernel calls reproduce the IRAF-style median.

5. IRAF parity tests and benchmarks

The repository’s benchmarks/IRAF workspace provides the IRAF comparison and throughput benchmark. Generated FITS inputs, IRAF outputs, manifests, CL scripts, and a local executable are ignored by Git. The parity test and benchmark normally create and remove temporary workspaces.

IRAF setup

Install the IRAF community distribution and make its ecl.e executable available at benchmarks/IRAF/ecl.e (a copy or symlink is fine). Set IMC_IRAF_ECL to use another executable and IMC_IRAF_ROOT to use another IRAF task tree. The benchmark also accepts --ecl and --iraf-root; an executable alone is insufficient because it needs the IRAF task tree. Without IMC_IRAF_ROOT, the benchmark discovers an adjacent ../iraf task tree when present. The parity test does not perform that lookup; set both environment variables for an external installation. The harness sets IRAFARCH=macos64; other IRAF builds are not covered by this configuration.

For manual inspection, first generate the inputs as shown below. Then run from the repository root with an explicit IRAF task tree and executable:

export IMC_IRAF_ROOT="/path/to/iraf"
export IMC_IRAF_ECL="$IMC_IRAF_ROOT/bin.macos64/ecl.e"
cd benchmarks/IRAF
env iraf="$IMC_IRAF_ROOT/" IRAFARCH=macos64 arch=.macos64 \
    "$IMC_IRAF_ECL" -f scripts/run_imcombine.cl

Do not add niterate to the generated local IRAF imcombine commands: the local imcombine.par does not define that task parameter. The general mapping from an IRAF-style iteration count to imcombiners remains maxiters, as shown in the table above.

Generated cases and parity test

Generate the inspectable local input matrix with:

uv run --extra bench python benchmarks/IRAF/scripts/make_inputs.py

It produces seven-frame, 128×128 uint16, int16, and float32 stacks by default. The unscaled cases use outtype=real, offsets=none, scale=none, zero=none, weight=none, project=no, and blank=0.0; the exact values for the rejectors are sigma=(2, 2), median centering, maxiters=5, nkeep=1, and revert_on_nkeep=True for sigma clipping; plus rdnoise=5, gain=2, and snoise=0 for CCD clipping; pclip=0.25; and n_minmax=(1, 1) for minmax.

The matrix covers sigclip, ccdclip, pclip, and minmax, with IRAF average/median (imcombiners mean/median) and these variants:

baseline: no pre-rejection threshold or input mask.
threshold: lthreshold/hthreshold of 990/1200 for uint16 and 90/220 for int16 and float32.
input_bpm: IRAF masktype="!BPM badvalue", maskvalue=1, mapped to a Boolean imcombiners mask.
snoise: unscaled CCDClip with snoise=0.05.
scaled: explicit per-image @file zero and scale values for CCDClip, PClip, and MinMax; CCDClip exercises IRAF’s sigscale branch.

Run the temporary-workspace parity suite with:

IMC_RUN_IRAF_PARITY=1 uv run --extra test --extra bench pytest tests/test_iraf_parity.py -q

It generates the inputs and CL script in a temporary directory, runs ecl.e, and compares the combined image, rejection masks, rejection counts, retained lower/upper bounds, and (except for the scaled case) the sigma diagnostic. Failures identify the case, IRAF parameters, imcombiners kwargs, differing entries, and example pixels. The suite contains a coverage check and 96 generated cases.

Benchmark

Run the benchmark with:

uv run --extra bench python benchmarks/IRAF/scripts/benchmark_iraf.py --repeats 3

It first validates output-only imcombiners results against IRAF, then reports batched all, sigclip, ccdclip, pclip, and minmax groups and, by default, one-case rows (--no-by-case disables the latter). Each IRAF workload repeats inside a single ecl.e -f process; its elapsed time is divided by the repeat count after subtracting the default dtype-specific 1×1 output-only imcombine baseline. Use --iraf-baseline startup or --iraf-baseline none to change that baseline.

The imcombiners timing reads the FITS inputs with fitsio, stacks them with NumPy, runs ndcombine(..., diagnostics=None), and writes the combined FITS output in a temporary directory. The resulting ratio is therefore a task-level comparison, not a standalone kernel microbenchmark. The script prints Markdown to standard output, including package, platform, and CPU information; use --output /tmp/imc-iraf-benchmark.md to retain a local report. Grouped rows are batched workloads rather than sums of the individual case rows. Published summaries belong on the 2-D Image Benchmarks page.