seiza (Python)

Project analysis

from pathlib import Path
import json
import seiza

report = seiza.analyze_project([Path("lights")], filter="Ha", samples=262144)
Path("analysis.json").write_text(json.dumps(report, indent=2), encoding="utf-8")
Path("analysis.html").write_text(seiza.analysis_html(report), encoding="utf-8")

analyze_project(inputs, *, samples=262144, filter=None, max_frames=None, metrics_only=False) accepts paths to files, nonrecursive directories or UTF-8 text lists; relative list entries resolve beside the list. max_frames limits sorted deduplicated paths before header filtering. It releases the GIL during analysis and returns a snake_case dictionary with schema_version, algorithm_version, options, frames, groups, timings, and warnings (schema/algorithm version 3). Invalid requests raise ValueError; individual inputs stay in frames with separate measurement_status and depth_status. Measurement failures use error; depth failures use depth_error. Measured frames contribute to native quality medians even if registration fails. metrics_only=True skips registration, normalization, sampling and replay; normalized noise and depth results remain unavailable. Groups count measured_frames, depth_frames and normalized_noise_frames separately. Use separate invocations per panel/pointing. No outcome is a quality verdict.

Mono, raw Bayer and RGB FITS/XISF frames are measured without producing a full stack. Per-filter groups contain registration/photometry diagnostics, common-footprint noise curves, known integration time, unknown-exposure counts and optional conditional trend estimates. analysis_html(report) returns a portable HTML dashboard with embedded data, filter controls and capture-time health plots. Neither function writes reports or modifies source images.

Noise measurements at even depths use balanced frame differences, which cancel shared systematic errors as well as scene structure. They do not reproduce production weighting/rejection or guarantee improvement with more exposure. See the complete method and schema.

balance_narrowband(image, shadow_protection=1, oiii_boost=1, sii_boost=1, max_samples=1000000) accepts interleaved display samples in [0, 1], with three physical planes H-alpha/OIII/SII, and returns (image, report) without changing the input. The report's channels list gives each physical channel's status, reason, samples, median, background, positiveDeviation and midtone. Curves protect the selected background; unreliable measurements are skipped. This nonlinear display operation is not photometric calibration. For HOO, duplicate OIII into the third plane and set has_sii=False; the third plane is ignored and the report includes only H-alpha and OIII.

screen_rgb_stars(nebula, original, starless) returns a new display RGB image with RGB stars screened onto the starless nebula. All inputs share one grid and lie in [0, 1]; original and starless RGB must share the same display transfer. Missing coverage masks whole pixels. Negative extraction noise is suppressed for display. The function does not run a star-separation tool.

Python bindings for seiza: star detection, WCS fitting, hinted/blind plate solving, satellite prediction, calibration, deconvolution, and batch/live image stacking for astrophotography, implemented in Rust.

pip install seiza

Binary wheels cover Linux (x86_64, aarch64), macOS (universal2), and Windows (x64); each is a single abi3 wheel for every CPython from 3.11 up. Type stubs are included, and computational image operations release the GIL. Some operations read input arrays in place while the GIL is released; others, including continuum fitting/subtraction and color composition, copy their inputs first. Do not mutate an array from another thread until the call returns.

Solve an image

The simplest way is to give solve_image the file. It reads the image, finds its stars, uses what the file says about itself (a WCS, the mount's pointing, the pixel scale, the capture time) and solves, blind when it must, as C hosts do with seiza_solve_image_json:

import seiza

solution = seiza.solve_image("M31.fits", catalog_dir="data", sip_order=3)
print(solution.source)          # "header-checked", "header-hint", "blind" or "header-unchecked"
ra, dec = solution.wcs.pixel_to_world(100.0, 200.0)

catalog_dir must directly contain the named star catalog .bin files; subdirectories are not searched. You can also pass a star catalog file. Header-checked and hinted solves need only that catalog. The blind .idx file is opened only when those stages cannot solve; an unavailable index then produces a specific fallback error, or preserves the existing header-unchecked result if the image contains a usable WCS.

FITS and XISF images have their stars measured on their own linear samples; PNG, JPEG and TIFF images are detected as decoded, with a photo's EXIF focal length as the first scale to try. To control each step instead, detect and solve separately:

import numpy as np
import seiza

# One-time: download the verified solver catalogs into the shared cache.
paths = seiza.fetch_catalogs()  # Tycho-2 solver + objects, Solar System, transients
catalog = seiza.StarCatalog.open(paths["stars-lite-tycho2.bin"])

# Detect stars in a 2D float32 (or uint8) luma array.
stars = seiza.detect(image_array)

# Hinted solve: approximate center and pixel scale. sip_order=3 also fits
# SIP distortion polynomials when enough matched stars support them.
solution = seiza.solve(
    stars, catalog, width, height,
    ra=150.1, dec=35.2, scale_arcsec_px=2.5, sip_order=3,
)
print(solution)                 # center, scale, matches, RMS
print(solution.rotation_deg, solution.flipped)
ra, dec = solution.wcs.pixel_to_world(100.0, 200.0)

open takes a file, a directory (the right catalog inside is picked — the deepest star catalog wins), or nothing at all. With no argument the standard places are searched: SEIZA_STAR_DATA / SEIZA_BLIND_INDEX, files next to the program, and the seiza setup directories (SEIZA_CATALOG_DIR). These are the same rules as the CLI's --data:

catalog = seiza.StarCatalog.open("data")   # directory
catalog = seiza.StarCatalog.open()         # after seiza setup

Star files hold positions at one epoch, 2025.5 for the published Gaia files. When a proper-motion sidecar sits beside one (stars-gaia.pm.bin beside stars-gaia.bin, written by seiza build-data gaia --motions), open reads it too, and catalog.has_proper_motions is true. Pass a solve the image's date as a Julian year, and it searches the stars where they were then:

solution = seiza.solve(stars, catalog, w, h, ra=..., dec=...,
                       scale_arcsec_px=..., epoch=2024.92)

Without a sidecar, or with epoch=None, the catalog is searched as written.

A solution is only as exact as the star positions it fits. When the array you detect on is stretched, or the sky has a gradient, measure the stars again on the image's linear samples before solving. refine_centroids fits each star's background as a plane over an annulus around it and iterates a Gaussian-weighted centroid, which is how the CLI and the C ABI place the stars of a FITS or XISF image. For a raw colour mosaic, pass it summed in 2 × 2 blocks, so every sample holds one red, one blue and two green pixels:

stars = seiza.detect(stretched_u8)
stars = seiza.refine_centroids(stars, raw_samples.astype(np.float32))
# A raw colour mosaic (height and width even), binned 2 x 2:
binned = raw.reshape(h // 2, 2, w // 2, 2).sum(axis=(1, 3), dtype=np.float32)
stars = seiza.refine_centroids(stars, binned, binning=2)

Stars can also be plain (x, y, flux) tuples from any other detector — the solver only needs positions and relative brightness:

solution = seiza.solve([(x1, y1, f1), (x2, y2, f2), ...], catalog, w, h,
                       ra=..., dec=..., scale_arcsec_px=...)

Measure stars and sensor tilt

The measurement detector is separate from the fast alignment detector above. It reports HFR, FWHM, SNR, flux, and optional Gaussian/Moffat PSF fits from a mono uint16 frame. The returned coordinates and radii are always in input pixels, including when detection binning is enabled:

measured = seiza.detect_measured_stars(
    image_u16,
    focal_length_mm=550.0,
    pixel_size_um=3.76,
    psf_type="moffat4",
)
cells, tilt = seiza.tilt_analysis(measured)
triangle = seiza.triangle_tilt_analysis(measured, angle_degrees=0)
print(len(measured.stars), measured.average_hfr, measured.average_fwhm)
print(tilt.tilt_percent, tilt.curvature_percent, tilt.worst_corner)
print(triangle.ready, triangle.tilt_percent, triangle.worst_sector)

The nine cells cover a 3×3 sensor grid. Fitted star theta and cell mean_theta values are ellipse major-axis orientations in radians over [0, π); theta_coherence describes how consistently stars in that cell share the direction. Tilt and curvature need measurements in all four corner cells, and curvature also needs the center cell.

tilt_analysis supplies the 3×3 measurements for a parallelogram diagram. triangle_tilt_analysis reuses the same detected stars and groups an inscribed circular annulus around three adjustment-screw axes. Its angle uses image coordinates: 0 points to the top and positive values turn clockwise. Sector IDs are 1-based, and their axes are the normalized input angle plus 0°, 120°, and 240°. Per-sector medians remain available for sparse data, but ready, tilt_percent, and best/worst sector enforce the native minimum of three stars per sector. overall_median_hfr is the median of every selected annular star, not the median of the sector medians.

detect_measured_stars_file(path, ...) measures a FITS or XISF file directly: its linear 16-bit luminance, with a Bayer mosaic debayered first, and with the focal length (FOCALLEN, FOCALLENGTH or FOCAL) and pixel size (XPIXSZ) read from its header when not given. Both calls take target_star_count, which retries with looser settings until that many stars are measured and otherwise returns the best pass.

For focusing tools, pixel_scale_arcsec(focal_length_mm, pixel_size_um) gives a rig's pixel scale, classify_pixel_scale(scale) the preset for it, and recommend_detection_binning(hfr) the detection_binning for a measured in-focus HFR.

Blind solve

No position hint, only a plausible scale range. Uses the prebuilt whole-sky pattern index and the deep Gaia catalog:

paths = seiza.fetch_catalogs(["stars-deep-gaia17.bin", "blind-gaia16.idx"])
catalog = seiza.StarCatalog.open(paths["stars-deep-gaia17.bin"])
index = seiza.BlindIndex.open(paths["blind-gaia16.idx"])
solution = seiza.solve_blind(stars, catalog, index, width, height,
                             min_scale_arcsec_px=0.5, max_scale_arcsec_px=15.0)

An image that already knows where it points can say so. read_solve_prior(path) reads a FITS or XISF file's own WCS (including a PixInsight solution), its RA/DEC or OBJCTRA/OBJCTDEC, and its pixel scale from PIXSCALE or pixel size over focal length, without decoding the pixels. solve_with_prior then checks that WCS, else searches around the position, else solves blind (given an index), and when nothing solves returns the file's own WCS with source == "header-unchecked":

prior = seiza.read_solve_prior("processed.xisf")
solution = seiza.solve_with_prior(stars, catalog, width, height, prior, index=index)
print(solution.source)   # header-checked, header-hint, blind or header-unchecked

A field that will not solve can search for a minute or more. Both solve and solve_blind take time_limit= in seconds, which raises TimeoutError, and cancel=, a callable that returns true to stop (asked about every 50 ms), which raises SolveError. Ctrl-C stops a solve either way, and an exception from cancel reaches you as itself:

solution = seiza.solve_blind(stars, catalog, index, width, height,
                             time_limit=30.0, cancel=stop_requested)

solution.false_match_probability is the chance that a wrong field would match as well, counting every candidate the search tried. A true field scores far below one in a billion; a blind solve above that is refused, and a hinted one above one in a thousand.

For faint fields, the optional stars-deep-gaia20.bin catalog reaches Gaia G≤20 (about 9 GB). It is intentionally not included in fetch_catalogs("all"), so request it explicitly with the same G≤16 blind index:

paths = seiza.fetch_catalogs(["stars-deep-gaia20.bin", "blind-gaia16.idx"])
catalog = seiza.StarCatalog.open(paths["stars-deep-gaia20.bin"])
index = seiza.BlindIndex.open(paths["blind-gaia16.idx"])

For phone and camera photos, read_photo_metadata reads the EXIF capture time, GPS position and focal length. Given the upright image's size (after PIL.ImageOps.exif_transpose), it also returns scale_ranges to try in order, the first from the focal length and the last a wide fallback:

from PIL import Image, ImageOps

image = ImageOps.exif_transpose(Image.open("phone.jpg"))
info = seiza.read_photo_metadata("phone.jpg", *image.size)
for min_scale, max_scale in info["scale_ranges"]:
    try:
        solution = seiza.solve_blind(stars, catalog, index, *image.size,
                                     min_scale_arcsec_px=min_scale,
                                     max_scale_arcsec_px=max_scale)
        break
    except seiza.SolveError:
        continue

FITS WCS output

Solutions convert directly to FITS WCS keywords (1-indexed CRPIX, TAN or TAN-SIP projection, CD matrix, and the complete A_p_q/B_p_q/AP_p_q/ BP_p_q set when distortion was fitted):

cards = solution.fits_header_cards()   # dict of keyword -> value
text = solution.fits_header_text()     # 80-column cards ending with END

The header text form is suitable for header-injection APIs — for example Siril's sirilpy scripting interface (set_image_header), which makes a seiza solve usable from a Siril Python script.

Objects in a solved field

field_objects lists the catalogued objects in a solved image: deep-sky objects and named stars, transients, and comets and asteroids, placed in its pixels with the outlines the catalog holds for them. It reads the catalogs fetch_catalogs installs (objects.bin, transients.bin, minor-bodies.bin); one that is missing gives no objects and an error that says why. read_capture_time reads when a FITS or XISF image (or a photo) was taken, which comets and asteroids need:

found = seiza.field_objects(
    solution.wcs, width, height,
    capture_time=seiza.read_capture_time("light.fits"),
)
for item in found.objects:
    print(item.name, item.kind, item.layer, item.x, item.y)
    for outline in item.outlines:
        print(" ", outline.role, outline.method, len(outline.contours))
print(found.object_catalog_error, found.minor_body_error)

Read pixel values from a file

read_image_window reads a rectangle of a FITS or XISF image's samples without decoding the rest of the file, so a viewer can show the values under a cursor on an image of any size. Only the rows of the window are read. The samples are in physical units: FITS samples have BZERO, BSCALE and BLANK applied (a blank sample is NaN), and XISF samples are as stored. Rows count in the file's stored order.

window = seiza.read_image_window("stack.fits", x=1200, y=840, width=16, height=16)
window.values.shape          # (planes, height, width), float32
window.x, window.y           # the window's corner, clipped to the image
value = window.values[0, 3, 5]   # plane 0 at (window.x + 5, window.y + 3)

A PixInsight XISF image normalized to 0..1 (bounds="0:1") reads as stored unless you pass normalized_full_scale, which puts it on that scale as stacking does: normalized_full_scale=65535.0 shows its values on the same 16-bit scale as camera frames. FITS images read the same either way.

master = seiza.read_image_window(
    "masterLight.xisf", x=1200, y=840, width=16, height=16,
    normalized_full_scale=65535.0,
)

A rectangle that runs off the image is clipped to it, and one wholly off it gives an empty window. A side may be at most 1024 pixels. A compressed, CIELab or complex XISF image cannot be read this way, because its samples are decoded as a whole; open it instead.

read_image_header reads a FITS or XISF file's header without decoding its pixels. headers maps each keyword to its value; an XISF file's PixInsight plate solution appears there as WCS cards. history and comments hold a FITS file's HISTORY and COMMENT text, where calibration and stacking programs record what they did.

header = seiza.read_image_header("light.fits")
header["headers"]["EXPTIME"]   # 300.0
header["history"]              # ["Calibrated with a master dark", ...]

Background extraction

Fit a compact background model to a C-contiguous mono (H, W) or RGB (H, W, 3) linear float32 array, inspect it, and then correct the image:

model = seiza.fit_background(stack, model="automatic", degree=2)
print(model.diagnostics)

corrected = model.correct(stack)                 # additive subtraction
illumination_corrected = model.correct(stack, mode="divide")
partial = model.correct(stack, strength=0.6)      # tune without refitting
background = model.render()                      # explicit full-size model

Fitting uses deterministic low-noise sample windows, robust sample rejection, and independent per-channel surfaces. Automatic mode chooses among polynomial degrees from held-out sample errors. Use model="radial_basis" with rbf_smoothing, or set allow_radial_basis=True in automatic mode, for an irregular thin-plate model. RBF is explicit because background samples can share real extended emission. model.correct() allocates only the corrected array; a full-size background exists only after render(). Pass a boolean (H, W) mask to exclude extended objects, dark clouds, registration borders, or source masks:

model = seiza.fit_background(stack, mask=structure_mask,
                             degree=1, samples_per_axis=12,
                             sample_radius=20)
for x, y, values, dispersion, weight, status in model.samples():
    print(x, y, values, status)

Fitted one channel at a time, the channels of a colour image can each pick a different automatic surface, and the mismatch shows as a colour cast toward the edges. select_shared_model chooses one surface from the held-out errors of every fit, and returns the fit_background keyword arguments to fit each channel again with, or None when the fits share no choice:

models = [seiza.fit_background(channel, model="automatic") for channel in channels]
shared = seiza.select_shared_model(models)       # the same automatic settings
if shared is not None:
    models = [seiza.fit_background(channel, **shared) for channel in channels]

The output remains linear and may retain negative or greater-than-one values. Background extraction is not display stretching or color calibration.

Light deconvolution

Apply the same conservative linear-image restoration as the Rust crate and CLI to a C-contiguous mono (H, W) or RGB (H, W, 3) float32 array:

restored = seiza.deconvolve(stack, psf_fwhm=3.1)

psf_fwhm is a measured unsaturated-star FWHM in pixels. The defaults use four damped Richardson-Lucy iterations and blend 35% of the estimate into the input. Pass masked=True for registered images whose missing border samples are NaN: the mask stays in the output and does not darken nearby data. Without it, non-finite samples raise seiza.EngineError. The returned array remains linear float32; no clipping or display stretch is applied. The operation releases the GIL. Inspect identical stretches for noise, rings, saturated-star failures, and field-dependent PSF mismatch before using a stronger iterations or amount.

Photometric colour calibration

calibrate_color fits per-channel gains for a linear RGB image from Gaia DR3 star colours: each star's instrumental colour is fitted against Gaia BP−RP, and the gains render a star of the white reference's colour (the Sun's by default) neutral, with the sky neutralized too.

import seiza
from astropy.io import fits

header = fits.getheader("stack.fits")
wcs = seiza.Wcs.from_header(header)
gaia = seiza.gaia_photometry_cone(56.75, 24.10, radius_deg=4.8)  # ESA archive; keep it
# Or offline, from the hosted catalog (about 460 MB, fetched once):
path = seiza.fetch_catalogs(["stars-gaia-photometry.bin"])["stars-gaia-photometry.bin"]
gaia = seiza.gaia_photometry_catalog_cone(path, 56.75, 24.10, radius_deg=4.8)  # file or directory
calibration = seiza.calibrate_color(rgb, wcs, gaia)               # rgb: (h, w, 3) float32
print(calibration.gains, calibration.red_fit, calibration.blue_fit)
calibrated = calibration.apply(rgb)

fit_bp_rp_range (default (0.0, 1.8)) sets the colour range of the stars the fits use, rejection_sigma (default 3) their outlier rejection, minimum_stars (default 20) the fewest stars a fit accepts, and calibrator_max_g the faintest Gaia G a calibrator may have (by default two magnitudes brighter than the faintest row given). calibration.measurements lists every star measured, with its position, g, bp_rp, R, G, B flux, and used, whether the red and the blue fit kept it.

calibrate_color_image and calibrate_and_stretch place the Gaia stars themselves, from rows or the offline catalog, and never fail because calibration cannot run: the outcome says why (no_solution, no_catalog, no_stars, too_few_stars, not_rgb, invalid_options, fit_failed) and the image is left uncalibrated. calibrate_and_stretch then display-stretches the image. Once calibration has run, an unlinked stage is linked instead, so the balance survives, and result.linked_stages lists those stages. For an image made from a solved frame, Wcs.transformed carries the frame's solution through the pixel mapping, counted from 0 with pixel centres on whole numbers; here a half turn on a 1300 × 1100 grid and a crop at (40, 60):

wcs = frame_wcs.transformed(((-1.0, 0.0), (0.0, -1.0)), 1299.0 - 40.0, 1099.0 - 60.0)
result = seiza.calibrate_and_stretch(
    rgb,
    wcs,
    catalog=path,  # or gaia=rows; None searches the installed catalogs
    stretch={"model": {"type": "auto-mtf", "target_median": 0.25, "shadows_clip": -2.8},
             "color_strategy": "linked", "max_analysis_samples": 200000},
)
print(result.calibration.status, result.calibration.reason)
display = result.image  # (h, w, 3) float32 in [0, 1]

A sensor clips each channel at its full scale, and a stack, scaling each frame's channels apart, leaves a clipped star core at a different level in each channel, which the colour gains then tint, often blue. A live stack records where its frames' raw pixels clipped: stacker.clipped_fraction() returns each pixel's clipped share of its observations, from 0 to 1, as a (h, w) array, or None when no frame clipped or the saturation level is unknown (it comes from SATURATE, SATLEVEL, or the frames' integer encoding). Pass it as clipped= to calibrate_and_stretch, or to seiza.neutralize_clipped(calibrated, clipped) after calibration.apply, and those pixels come out neutral: each pixel's channels rise toward its brightest, fully from a share of one half, and the map is feathered over a few pixels so a white core fades into its star's halo. A drizzled stack takes the live stack's map as it is; each of its pixels reads the map pixel it lies in.

Image stacking

The wheel includes the same linear calibration, registration, normalization, and online rejection engine as the Rust crate and CLI. Batch stacking accepts FITS paths and writes an unstretched linear float32 FITS result:

options = seiza.StackOptions(
    normalization="local-background",
    local_tile_size=256,
    demosaic="vng",
    maximum_drift_pixels=256.0,
    maximum_drift_fraction=0.15,
)
result = seiza.stack_fits(
    sorted(light_paths),
    "stack.fits",
    options=options,
    bias="master-bias.fits",
    dark="master-dark.fits",
    flat="master-flat.fits",
)
for frame in result.frames:
    print(frame.source, frame.accepted, frame.reason, frame.registration_rms_pixels)

normalization is "none", "global" (the default), "local", or "local-background", which keeps one gain per channel and matches each tile's background to the reference's to reduce seams. Offsets use the median of differences at corresponding sky positions, after applying the gain, so shared nebular structure cancels before the background fit. On an expanding mosaic canvas, which keeps the first value each sample received, the live pass fits one background offset per frame instead, so the stack does not take on the sky of whichever frames filled the canvas first. Mosaic reintegration builds this reference from the best exposures covering each part of the sky, up to twenty per sample. It preserves their recorded spatial corrections and fits remaining additive differences, keeping the gains and shared sky structure. demosaic chooses how one-shot-colour frames are demosaiced: "vng" (the default) keeps each star's colour across its profile, "mhc" gives slightly sharper stars but rings around small ones, and "bilinear" is fastest. registration_model is "similarity" (the default here), "affine", or "quadratic", which fits a polynomial to follow lens distortion, for example across a meridian flip; the CLI defaults to "quadratic". weighting is "equal" (the default) or "inverse-noise-variance", which weights each frame by the reference frame's noise over its own, squared. The CLI's Bayer drizzle option is not yet exposed here.

Under noise weighting, each accepted frame's noise and weight list one value per channel, and LiveStacker.reference_noise gives the reference frame's noise, the unit of the weights; under equal weighting they are empty. seiza.frame_noise(image) measures any linear array the same way. Each accepted frame's mapping, and LiveStacker.reference_mapping for the reference, is a dict of its registration and normalization, to keep beside the frame. Superpixel mappings include channel_sites, three [x, y] offsets in source pixels in RGB order. Saved contexts retain these offsets, so crop extraction from equally prepared source frames keeps colours aligned.

disposition.resample(image, interpolation="lanczos3") resamples another image of the frame's size onto the stack's grid by that frame's registration alone, without its normalization: its dark, a mask or a defect map. It resamples as the stack did, warp and channel sites included, and leaves NaN where the image does not reach. A rejected frame has no mapping to resample by, and raises ValueError.

LiveStacker.open_prepared_fits(path) opens a stack from a frame already calibrated and debayered, keeping its FITS header cards, so you can correct a frame between calibration and registration; push later frames as arrays.

Without a dark, stack_fits learns the sensor's hot and dead pixels from up to 16 of the lights and repairs them in every light before debayering and registration, as the CLI's --defects auto does: a defect stays on its pixel while dithering or drift moves the sky, so a pixel counts only when it stands out in most of the lights, and never when its neighbours rise as a star's do. defects="on" does so with a dark as well, for the defects its master leaves, and defects="off" never. seiza.detect_defects(paths, bias=..., flat=...) gives the map on its own, a DefectMap with width, height, frames (the lights examined) and pixels, raw-frame indices y * width + x. A steady hot pixel's excess over its neighbours, learned at exposure_seconds, is subtracted from each light as a dark would subtract it, so a star the dither carries onto it keeps its core; dead and wandering pixels take the median of their neighbours. Pass the map as StackOptions(defects=...) to any stack, or build one from a known list with DefectMap(width, height, pixels). Its sigma (5), persistence (0.6, the fraction of lights) and maximum_frames (16) are the CLI's.

The first light is the reference unless reference="auto", which scores every light and takes the best, as the CLI's --reference auto does. The result's reference_light names the light used. seiza.choose_reference(paths) gives the same choice on its own: its reference_path, a ranking of the scored lights' indices, best first, and each light's scores.

To stack each filter onto the same grid, pass every run the same reference_image, as the CLI's --reference-image does. It sets only the alignment grid and WCS; the reference light still sets normalization and noise weighting, and every light listed is stacked, that one included. With reference="auto", if the best light does not align to the reference image, the next-best lights are tried in turn, up to five in all. The reference image is not calibrated and adds no pixels unless it is also listed as a light. It must have the lights' channel count:

for name in ["red", "green", "blue"]:
    seiza.stack_fits(lights[name], f"{name}.fits", reference_image="common.fits")

LiveStacker.open_fits_with_reference_image(reference_image, channel_reference) does the same for a live stack. It starts with no frames, so push the channel reference like every other light; reference_integrated is then false. LiveStacker.open_fits_with_reference_image_candidates(reference_image, channel_references) takes several lights, best first, and uses the first that aligns. It returns the stack and a selection naming the light used and why each earlier one was skipped:

choice = seiza.choose_reference(red_lights)
best = [red_lights[index] for index in choice.ranking[:5]]
stacker, selection = seiza.LiveStacker.open_fits_with_reference_image_candidates(
    "common.fits", best
)
for index, path, reason in selection.skipped:
    print(f"{path} did not align: {reason}")
stacker.push_fits(selection.path)

For separate filter mosaics on one grid, first plan an expanded stack covering the desired field, then use its open_fits_on_canvas method:

options = seiza.StackOptions(normalization="global")
ha = seiza.LiveStacker.open_fits_expanding(ha_paths, options=options, flat=ha_flat)
oiii = ha.open_fits_on_canvas(oiii_paths, options=options, flat=oiii_flat)
sii = ha.open_fits_on_canvas(sii_paths, options=options, flat=sii_flat)
for stack, paths in [(ha, ha_paths), (oiii, oiii_paths), (sii, sii_paths)]:
    for path in paths:
        stack.push_fits(path)

Each stack starts empty. Only its own lights contribute pixels, calibration, normalization and noise weighting. The canvas's dimensions, panel geometry, lens correction and WCS (when available) are shared. A stack reopened with open_context can also supply the canvas. Saving a context requires accepting at least one light first.

open_fits_on_canvas(paths, *, options=None, bias=None, dark=None, flat=None, dark_exposure_seconds=None, cancel=None) never grows the source canvas: samples outside it are cropped and missing channel coverage stays empty. Plan the full desired footprint first, and use the same bin factor for every channel. The first supplied light anchors that channel's brightness; lights without photometric overlap are retried after connecting panels join. Failed joins are rejected when pushed. Cancellation and input/calibration errors abort construction without changing the source stack.

LiveStacker.open_fits_expanding(paths) grows the canvas past the anchor to hold every connected, overlapping light, as the CLI's --reference expand does: mosaic panels, including ones that connect only through other panels, or a dithered stack whose edges should be kept. The first path anchors the WCS and photometric scale, so put choose_reference's choice first to anchor on the clearest light. Lights may differ in rotation by any amount. A light that cannot join, such as one of another target, comes back rejected with the reason; StackError is raised only when nothing joins the anchor. The stack starts empty; push every light, the anchor included:

from pathlib import Path

anchor = seiza.choose_reference(panel_lights).reference_path
order = [anchor] + [path for path in panel_lights if Path(path) != anchor]
options = seiza.StackOptions(
    normalization="local-background",
    registration_model="quadratic",  # warps frames across a meridian flip
)
stacker = seiza.LiveStacker.open_fits_expanding(
    order, options=options, dark="master-dark.fits", flat="master-flat.fits"
)
dispositions, report = stacker.push_fits_pipelined(order)
stacker.finish("mosaic.fits")

max_canvas_megasamples (256 by default) bounds the canvas's pixels times channels; a run peaks at about 40 bytes per sample. When the panels' overlaps show one lens distortion bending them all, planning undoes it and stacker.expanded_lens describes it as a dict (its axis, coefficients, largest shift in pixels, and the overlaps' residual with and without it); otherwise it is None. stacker.expansion_planning_notes lists what planning noted: first, when lights did not join the canvas, how many, and the most that those matching a panel overlapped it against minimum_overlap; then why it skipped fitting one lens to every panel, or why it kept its first plan when the lens-corrected one failed. Unless options set minimum_overlap, a light joins the canvas sharing 0.05 of its footprint with it, as with the CLI's --reference expand; a fixed reference needs 0.60. A mosaic planned for 10% overlap lands nearer 9%, and a minimum of 0.1 can leave every panel but the anchor's off the canvas. cancel, a callable, is asked before each light is planned, and Ctrl-C also stops planning. Drizzle is not supported on an expanded canvas. See the CLI documentation for how panels are matched and what the result looks like on real data.

stacker.pipeline_memory gives what preparing frames needs, frame by frame, for planning memory: worker_bytes for each worker, integration_bytes for the frame being integrated, and sequential_bytes for one frame at a time. A frame joining an expanded canvas is resampled onto only the part of the canvas it covers, so it is costed at the largest frame planning read, at any rotation and scale the acceptance admits, rather than at the canvas's size; a derived workers count follows the same sizes. On a four-panel mosaic of 88 megapixels from 26-megapixel frames, each worker is costed at 2.4 GiB rather than the canvas's 6.5 GiB.

Lens-corrected mosaics preserve each panel's sky projection even with the default registration_model="similarity", including new exposures added after planning or reopening a saved context. Older saved lens-corrected contexts without panel projections can still be replayed; replan from the source lights before adding new exposures.

For live integration, construct from a FITS path or a C-contiguous mono/HWC RGB NumPy float32 array. An array-based stacker accepts only already-linear, calibrated arrays through push(). A path-based stacker can use push_fits() for its configured calibration path or push() for caller-prepared arrays. A stacker keeps that input mode after checkpointing. Both methods return a typed admission decision, and a rejected frame never mutates the accumulator:

stacker = seiza.LiveStacker.from_array(reference, options=options)
for frame in incoming_arrays:
    disposition = stacker.push(frame)
    if not disposition.accepted:
        print(disposition.reason)

preview_state = stacker.snapshot()  # immutable copy
linear_mean = preview_state.image
coverage = preview_state.coverage
final = stacker.finish("stack.fits")  # consumes the live accumulator

To save the stack so far and keep going, call stacker.write_fits(path), which copies only the mean, or write_fits(path) on any StackSnapshot, including one from reintegrate(). Both write linear float32 FITS (XISF for a .xisf path) with the reference frame's WCS, and refuse a path that is one of the stack's inputs or masters.

stacker.state() describes the stack: its size and counts, a configuration_fingerprint of its options, masters and input mode (to check that a resumed stack is the one you meant), its input_mode, the input_paths it has used, its reference_frame (role, master flag, signature and calibration state), its reference_image_path, and reintegration_unavailable, which says why reintegrate() would fail, or is None. stacker.set_calibration_fits(bias=..., dark=..., flat=...) replaces the masters for frames pushed from then on; frames already stacked keep the calibration they had. Every file is read and checked against the reference before anything changes, and passing none clears the masters:

stacker.set_calibration_fits(dark="master-dark-300s.fits", flat="master-flat.fits")
print(stacker.state().configuration_fingerprint)
stacker.write_fits("so-far.fits")

When the paths are all known — a finished session rather than a live one — push_fits_pipelined() prepares several frames at once while integrating in the order given, so the result is identical to pushing them one at a time. It measured 2.0x faster on local storage and 3.0x with a 300ms read latency:

dispositions, report = stacker.push_fits_pipelined(paths)
print(report)  # PipelineReport(integrated=..., rejected=..., failed=...)

A path that cannot be read, or that repeats one already stacked, comes back in place with accepted false and a reason rather than raising, so one bad path in a night's listing does not lose the rest — check report.failed rather than reading a clean return as success. Pass workers= when the frames arrive over a network, since the library cannot tell a network mount from a local disk, and normalized_full_scale=65535.0 when the set mixes PixInsight XISF frames with 16-bit camera data.

Before pushing a light under masters you chose yourself, ask whether they accept it. validate_calibration() raises StackError naming the reading that failed and the tolerance in force; compatible_calibration() answers the same question as the masters it would keep and why it would drop the rest. Pass the MatchTolerances you matched the masters under, so both judge with your gate:

light = seiza.FrameSignature(width=6248, height=4176, channels=1, rotation_deg=113.8)
wide = seiza.MatchTolerances(rotation_deg=180.0)
stacker.validate_calibration(light, wide)        # None, or raises StackError
kept, dropped = stacker.compatible_calibration(light, wide)

For flats taken at a different gain or offset, opt in with MatchTolerances(allow_flat_gain_offset_mismatch=True) and pass it as StackOptions(match_tolerances=...). Use this only after verifying that the camera's normalized flat response transfers across those settings. The loaded master flat must record DARKSUB=T; normalization alone does not qualify. Build it with flat-darks matching the flats' gain, offset, readout mode and exposure. Light darks still need to match the lights. The option also permits selecting different-gain flats in plan_calibration, but never mixes different gains within one master or relaxes dark-flat selection. Saved contexts and reintegration retain the policy.

StackOptions(maximum_registration_rms=2.0) sets an absolute registration limit in original sensor pixels. There is no running-median RMS gate. Plate-solve RMS measures a different fit and does not establish that a similarity transform aligns mosaic panels. registration_model="quadratic" can model curvature a similarity transform cannot; inspect alignment before choosing it. interpolation="lanczos3" resamples each frame through Lanczos-3, clamped against ringing at high-contrast edges: sharper than the default "bilinear", and the CLI's default. Rejections include the actual reason. Fraction thresholds accept equality, and their messages preserve the precision of the measured fraction.

Passing tolerances replaces the stack's own whole: a field you leave unset takes the built-in default, not the stack's value. Both calls raise StackError for a tolerance that is NaN, infinite, negative, or a rotation above 180 degrees, which already accepts any angle.

The stack judges every light it is given under its own gate, so give it the same one, or the light is refused when you push it. StackOptions(match_tolerances=...) sets the gate a stack opens with, which also judges the reference. set_match_tolerances() changes it for the lights pushed after, and match_tolerances() reads it back. Each admitted frame keeps the gate it was admitted under, and reintegrate() judges it by that gate, so you can widen the gate for one light and narrow it again after:

stacker.set_match_tolerances(wide)
stacker.push_fits("light-113.8deg.fits")
stacker.set_match_tolerances()                   # the defaults again

A rotator's angle can be wrong for the question a flat asks: after a re-home it can read the same angle for another position. When you have measured which flats a light's pixels follow, by plate solving it or fitting its sky against the candidate flats, give that angle as light_rotation, and the gate compares it in place of the light's recorded angle. "unknown" compares the light at no angle, so any flat angle passes; "recorded" is the default. The flats keep their recorded angles, and plan_calibration reads it the same way:

stacker.set_match_tolerances(seiza.MatchTolerances(light_rotation=3.5))
stacker.match_tolerances().light_rotation        # 3.5

A stacker built from FITS or XISF paths can integrate its admitted frames again from their files, with leave-one-out rejection that removes trails the online pass kept, and can drizzle that integration onto a finer grid as WBPP's DrizzleIntegration does. Bayer frames drizzle each photosite into its own colour:

replayed = stacker.reintegrate()            # StackSnapshot; the live stack is unchanged
integrated, drizzled = stacker.reintegrate_drizzled(scale=2)
drizzled.write_fits("drizzle-2x.fits")      # reference WCS scaled to the grid
weights = drizzled.weight                   # drop area times frame weight

Each of those reads every frame from its file and prepares it again. Calling stacker.retain_frames_for_reintegration(scratch_directory=...) before pushing frames keeps each admitted frame's registered image on disk as it is stacked, about four bytes per pixel per channel for each frame, so a reintegration reads that instead and gives the same result sooner.

reintegrate(frames=paths) integrates only the admitted frames read from those files, in that order. Each is normalized and weighted as in a reintegration of every frame, against the same background reference, so stacks of disjoint subsets share one grid and one normalization and differ only by their frames. With retained frames, one registration serves any number of subsets:

stacker.retain_frames_for_reintegration(scratch_directory="/nvme/scratch")
...                                         # push every frame of the night
half_a = stacker.reintegrate(frames=paths[0::2])
half_b = stacker.reintegrate(frames=paths[1::2])
noise = half_a.image - half_b.image         # independent noise, same sky

A path the stack did not admit, an empty list or a frame named twice raises StackError.

For large mosaics, pass an existing writable directory on a disk with room as scratch_directory to reintegration. An unusable explicit directory is an error. With the implicit system temporary directory, cache failure may fall back to reading whole frames and require substantially more RAM. Mosaic scratch images and their background sample files store only the registered footprints and intersecting sample tiles.

StackOptions(bin_factor=2) (or 3, or 4) area-mean bins each light after it is calibrated, repaired and debayered, before registration, so a large sensor or a mosaic stacks quickly and in a quarter (or a sixteenth) of the memory. Pixel options keep their meaning in the sensor's pixels, and registration_rms_pixels reports in them; the live view, snapshots and each frame's mapping (which then carries "bin_factor") are on the binned grid, and measure_depth() says what the stack was binned by, since binning lowers each pixel's noise. Bayer drizzle needs single photosites, so it is refused with a bin factor above 1. With an even bin factor, demosaic="superpixel" makes each 2x2 cell of a colour sensor one pixel, about ten times faster than VNG; registration reads each colour from where its photosite lies, so the stack's colours line up. When the session is done, reintegrate the same frames at full resolution without registering them again:

options = seiza.StackOptions(bin_factor=2, normalization="local-background")
stacker = seiza.LiveStacker("light-001.fits", options=options)
...                                         # push the rest as usual
full, refined, rms = stacker.reintegrate_full_resolution()
full.write_fits("stack-full.fits")          # WCS carried to the sensor's grid

Each frame's transform is carried exactly from the binned grid, and with refine=True (the default) a single field refits it against stars detected at full resolution, as registration would find them there; a mosaic refits unwarped frames against the canvas the other frames build at full resolution. Frames with a panel warp keep it; when all frames have warps, the unnecessary refit and its input reads are skipped. refined counts the frames refitted and rms is the median registration RMS in the sensor's pixels, or None when no frame's is known. A superpixel stack is demosaiced with VNG at full resolution. For a quicker look, reintegrate_drizzled(scale=2) on a stack binned by 2 drops each frame's sensor pixels onto the sensor's grid, with rejection decided at the binned resolution.

Frames calibrated, registered and normalized somewhere else, one size on one pixel grid with uncovered samples as NaN, take the same rejection through integrate_registered_frames. It reads each file once, keeps its samples in a scratch file inside scratch_directory (the system temporary directory by default; four bytes per sample, removed when it returns), runs the three passes over those files in bands of rows, and registers nothing again; frame_weights gives one weight per channel for each frame:

stack, frames = seiza.integrate_registered_frames(registered_paths, sigma_high=3.0)
stack.write_fits("stack.fits")                   # with the first frame's WCS
for frame in frames:
    print(frame.path, frame.finite_samples, frame.integrated_samples)

For registered mosaic panels, pass mosaic_panels=[0, 0, 1, 1] when the first two inputs cover one panel and the next two cover another. Seiza matches additive background planes across overlaps, then tapers remaining smooth differences at panel edges. Inputs must share a pixel grid and photometric scale; missing pixels stay NaN. The option accepts at most 256 distinct panel IDs and leaves the existing rejection and frame weights in use. Matching compares paired positions before clipping, including partial overlap cells and independent channel coverage. See mosaic background matching.

Where the hardware has a GPU, a build with seiza-stacking's gpu feature (the default) runs stacking's heavy passes there, in single precision, so results differ from the CPU's by rounding. A GpuOptions says which passes may: give it as gpu to StackOptions, which governs the stack and its reintegration, to detect_defects and to integrate_registered_frames. GpuOptions(enabled=False) keeps every pass on the CPU; defects, debayer, detection, resampling, accumulator and integration keep one each. concurrency (3) is the frame calls the GPU takes at once. memory_bytes caps the device memory the GPU's work is planned in, by default what the device has: a live stack's accumulator stays on the CPU when it, the frames prepared for it and the frame calls' buffers would need more than three quarters of it. trace=True prints each call's timings to stderr:

options = seiza.StackOptions(gpu=seiza.GpuOptions(debayer=False))
stack, frames = seiza.integrate_registered_frames(
    registered_paths, gpu=seiza.GpuOptions(enabled=False)
)

A context keeps the stack's GPU options, and its configuration fingerprint leaves them out. LiveStacker.open_context(path, gpu=...) reopens a stack under other options, for example when a user turns the GPU off, and the fingerprint stays the same:

stacker = seiza.LiveStacker.open_context(
    "m31.seiza-stack", gpu=seiza.GpuOptions(enabled=False)
)

Checkpointing is non-consuming. Reopening preserves the original registration reference, calibration and options, online rejection statistics, coverage, and the FITS/XISF source ledger:

stacker.save_context("m31.seiza-stack")

stacker = seiza.LiveStacker.open_context("m31.seiza-stack")
decision = stacker.push(next_array)
stacker.save_context("m31.seiza-stack")

The context file is versioned, compressed, checksummed, and atomically replaced. It is processing state for resumption; finish("stack.fits") remains the interoperable final image output.

Frames taken after a German-equatorial-mount meridian flip are handled by default. maximum_rotation_degrees limits deviation from either the reference orientation or its 180-degree counterpart; frame diagnostics still report the full fitted rotation.

Snapshot array properties are copies, so Python cannot mutate live Rust state. All expensive FITS, calibration, registration, and integration work releases the GIL.

Color from mono stacks

Aligned mono float32 arrays can be combined without writing intermediate files. Outputs have shape (height, width, 3):

rgb = seiza.combine_rgb(red, green, blue)
lrgb = seiza.combine_lrgb(luminance, red, green, blue,
                          luminance_weight=1.0)
super_lrgb = seiza.combine_lrgb(luminance, red, green, blue,
                                luminance_mode="super")
super_rgb = seiza.combine_rgb(red, green, blue, luminance_mode="super")

sho = seiza.combine_narrowband(ha, oiii, sii, palette="sho")
hoo = seiza.combine_narrowband(ha, oiii, palette="hoo")
foraxx = seiza.combine_narrowband(ha, oiii, sii, palette="foraxx-sho")
custom = seiza.combine_narrowband_matrix(
    ha, oiii, sii,
    matrix={"red": {"sii": 0.7, "ha": 0.3}, "green": {"ha": 1.0},
            "blue": {"oiii": 1.0}},
)

combine_narrowband_matrix mixes the filters through your own weights: each output channel maps "sii", "ha" and "oiii" to a weight, and a weight left out is zero. A matrix that weights SII needs sii.

Each combine_* also writes its result when given output=: 32-bit float FITS, or XISF for a .xisf name, labelled in SEIZACLR and with SEIZATRF saying whether the samples are linear or display-referred. Pass the first channel's FITS header (an astropy Header or a dict) as header= and the output keeps its WCS, moved by any crop:

from astropy.io import fits

red, red_header = fits.getdata("red.fits", header=True)
seiza.combine_rgb(red.astype(np.float32), green, blue,
                  output="rgb.fits", header=red_header)

Pass crop="bounds" or crop="inscribed" to trim the blank edges that registering one channel onto another leaves behind. bounds keeps the box every channel covers; inscribed keeps the largest rectangle they all cover in full, so nothing stays NaN. seiza.crop_report measures the same thing without composing, and flags any channel whose coverage sits more than report["off_center_limit_pixels"] from the others:

report = seiza.crop_report({"red": red, "green": green, "blue": blue})
x, y, width, height = report["region"]
stray = [c["name"] for c in report["channels"] if c["off_center"]]

The default percentile normalization is a quick-look channel match. Pass normalization="none" for already matched inputs. Foraxx inputs must also already lie in [0, 1] in that mode; keep percentile normalization for sensor-unit arrays. RGB, LRGB, additive super-LRGB (L + R + G + B), synthetic super-RGB (R + G + B), the six direct S/H/O permutations, and HOO are linear-light. Super-luminance output can exceed one. Foraxx-SHO/HOO use a stretched working copy as required by the published dynamic formula, so those returned arrays are display-referred. Composition releases the GIL.

Parameterized display stretching

seiza.stretch applies the shared seiza-stretch model to mono (H, W) or RGB (H, W, 3) float32 arrays and returns display-referred float32 without eight-bit quantization:

preview = seiza.stretch(linear, model="percentile-asinh",
                        black_percentile=0.01,
                        white_percentile=0.995, strength=10)
preview = seiza.stretch(linear_rgb, model="auto-mtf",
                        target_median=0.2, shadows_clip=-2.8,
                        color_strategy="luminance-preserving")
preview = seiza.stretch(linear, model="ghs", stretch_factor=4,
                        local_intensity=-1, symmetry_point=0.35,
                        protect_shadows=0.1, protect_highlights=0.8)
preview = seiza.stretch(linear, model="median-anchor", target_median=0.25,
                        black_point_sigma=2.7, hdr_amount=0.3)
preview = seiza.stretch(linear_rgb, model="layered", target_background=0.15,
                        contrast_recovery=0.6, dynamic_range_compression=0.5,
                        saturation_amount=0.4)

Available models are identity, linear, asinh, percentile-asinh, mtf, manual ghs, auto-mtf, median-anchor, and layered; color strategies are linked, unlinked, and luminance-preserving. median-anchor places the background median exactly on target_median (default 0.25; auto-mtf defaults to 0.2), with the black point black_point_sigma robust sigmas below the median (None uses the image minimum), then optional curves_boost, hdr_amount and hdr_knee highlight compression, and normalize. layered stretches luminance with the same curve, target_background (default 0.15) and a black point set by aggressiveness (0 to 1, default 0.5), and treats large and fine structure apart: scale_separation (8 to 1024 pixels, default 256) splits them, contrast_recovery (0 to 1) keeps detail on bright structure, dynamic_range_compression (0 to 1) shrinks star cores, and saturation_amount, saturation_boost (both 0 to 1) and lightness_mask lift colour. It reads the array's height and width, and carries colour by luminance, so it takes linked or luminance-preserving. Analysis and application release the GIL.

Calibration masters use the same bounded-memory two-pass builder:

bias = seiza.build_bias(bias_paths, "master-bias.fits")
dark = seiza.build_dark(dark_paths, "master-dark.fits",
                        bias="master-bias.fits")
flat = seiza.build_flat(flat_paths, "master-flat.fits",
                        bias="master-bias.fits",
                        dark_flat="master-dark-flat.fits")

Each builder reads every input twice, so a master over dozens of frames runs for minutes. Pass cancel= a predicate to stop one early — it is called once per input and raises StackError when it returns true. Ctrl-C is honoured at the same points, without a predicate.

stop = threading.Event()
dark = seiza.build_dark(dark_paths, "master-dark.fits", cancel=stop.is_set)

build_flat(..., flat_star_masking=seiza.FlatStarMasking()) leaves stars and broad saturated patches out of each sky flat before it is normalized, and scratch_directory= keeps the build's temporary pixels in a directory you choose. build_dark(..., dark_level_screening=1.0) sets aside darks whose level sits more than that many times the quietest dark's pixel noise above it, the mark of stray light. The result's masked_samples, input_masked_samples and flat_star_masking say what masking left out, and rejection_method which rejection the build used.

A build sets aside an input whose header contradicts the first one, such as a flat from another rotator angle, rather than failing. result.skipped_inputs lists each as (path, reason). result.input_paths lists the inputs used, and result.input_statistics their (accepted, rejected) sample counts, in the same order. The build checks that these account for every input before it writes the master.

build_flat(..., defect_suppression=True) replaces hot and cold pixels in the flat master with their neighbours' median. A defective pixel repeats in every flat, so rejection across frames keeps it. defect_sigma_low and defect_sigma_high (16 by default) set how many sigma from its neighbourhood a pixel must sit; the default leaves star cores alone. result.defect_pixels_replaced counts the pixels it replaced.

Choosing calibration frames

seiza.probe_frame(path) reads a FITS or XISF header, without its pixels, into the frame's role ("bias", "dark", "dark-flat", "flat", "light" or "unknown"), whether it is_master, its FrameSignature and its calibration_state. seiza.plan_calibration applies the whole matching policy to those probes, as the desktop apps do. It keeps the candidates that suit every light, orders them by capture time, nearest the reference light first, and takes the first set of minimum frames that agree with each other. It says why it left out each of the rest:

light = seiza.probe_frame("light-001.fits")
darks = [seiza.probe_frame(path) for path in dark_paths]
plan = seiza.plan_calibration("dark", light, darks, minimum=5)
if plan.ready:
    seiza.build_dark(plan.selected_paths, "master-dark.fits")
for path, reason in plan.excluded:
    print(path, reason)  # for example "exposure-mismatch"

kind is "bias", "dark", "dark-flat" or "flat". Dark flats are planned against the flat they will calibrate rather than a light. Pass references= to choose one set for several lights, tolerances= a MatchTolerances, and bias_available=True once you have a bias master, so that darks of another exposure can be scaled. A seiza.CalibrationRecord(path, signature, role=...) stands in for a probe when your records come from elsewhere. A bad request raises ValueError.

Screening calibration frames

Before a master is built, judge its frames as a set. measure_dark and measure_flat read one raw FITS or XISF frame each and return a small measurement; keep it (to_json(), DarkMeasurement.from_json()) and judge it again whenever its set changes. screen_bias, screen_darks and screen_flats give each frame a FrameVerdict: a grade ("ok", "suspect" or "poor") and the findings behind it, each with its code, value and the limit it reached. Most checks compare a frame with its peers and need three or more frames. The verdicts say what the pixels show, never why:

darks = [seiza.measure_dark(path) for path in dark_paths]
for path, verdict in zip(dark_paths, seiza.screen_darks(darks)):
    print(path, verdict.grade, [finding.code for finding in verdict.findings])

flats = [seiza.measure_flat(path) for path in flat_paths]   # raw, before calibration
report = seiza.screen_flats(flats, clipped=(0.001, 0.01))   # (suspect, poor)
print(report.master.grade, report.changing_features)

A flat is measured against full_scale (65535 by default), and is read as a raw colour frame when its header records a Bayer pattern (cfa= decides otherwise). Each limit is a (suspect, poor) pair.

Residual flat responses

A flat cannot correct dust that moved after it was taken. Given crops of several lights cut at the same place on the detector, build_residual_flat keeps only the light most of them lose, and returns that response with its evidence. Show first that the dip stays put on the detector while the sky moves; this does not decide what it is:

response, diagnostics = seiza.build_residual_flat(crops)   # one crop per light
seiza.validate_residual_flat_response(stored_response)      # raises ValueError
corrected = seiza.apply_residual_flat_response(response, light, x, y)

The response is 1 where nothing is corrected and never gains more than maximum_gain (1.2). apply_residual_flat_response returns a corrected copy with the response's top-left corner at pixel (x, y).

Flat sessions and their dust

A flat set reduces to one response: its flats binned and combined. Two responses of one light path compare by their ratio, and the change between them carries a flat across, as when dust moved between two sessions. The same calls as the C ABI's seiza_flat_* and seiza_dust_*; planes are 2-D float32 arrays, and results are dicts of the core's own JSON:

factor = seiza.flat_bin_factor(width, height)
planes = [seiza.flat_bin_frame(flat - pedestal, factor) for flat in flats]
response = seiza.flat_build_response(planes, factor)        # {"plane", "level", ...}
comparison = seiza.flat_compare(earlier["plane"], later["plane"], factor)
change = seiza.flat_change_map(earlier["plane"], later["plane"])
seiza.flat_apply_change_map(light, change, factor, divide=True)   # in place

session = seiza.dust_find_session(raw_planes, factor, pedestal)   # dust, counted per flat
seiza.dust_frame(raw_planes[0], pedestal, session)                # one flat against it
seiza.dust_find_changes(earlier["plane"], later["plane"], factor) # what came, went, moved

flat_lights_background, flat_fit_change and flat_choose_response weigh a few lights' own sky against a change or a set of candidate responses; flat_encode_plane and flat_decode_plane keep planes in a cache. Options are a dict of DustOptions fields, such as {"edge_sigma": 3.0}.

FlatGroups holds PSF Guard's measured grouping policy: response families, groups, changes, angle epochs, uncertain light placements and resolved personal choices. Build it once from the stored rows, a mapping of responses, comparisons, evidence, imaging (or imaging_lights with profile_rigs), limits, group_overrides, night_boundaries, epoch_samples, family_evidence and light_choices; omitted collections are empty. Its properties are the book (families, boundaries, events, overrides, outliers, spanning, angle_terms, epochs), and pairs_to_compare, sessions_in_windows and sessions_needing_evidence are the measurement jobs and evidence keys for later pixel passes. place_light(light) says where one light stands: its automatic family among the compatible_families it names, the pixel question that would settle a close call, its resolved choices, and its placement against each change. Records use snake_case; times are Unix seconds, angles degrees, levels ADU, dust/broad fractions of background. Compatible families must already pass the host's sensor/optics gates. Automatic results and personal choices stay distinct. This stage does not grade frames or build masters. See the schema and executable fixture.

groups = seiza.FlatGroups(rows)
groups.families[0]["groups"]
placed = groups.place_light({"file_name": "L_0001.fits", "rig_uuid": rig,
                             "filter": "Ha", "at": captured_at, "angle": 12.0,
                             "compatible_families": [groups.families[0]["key"]],
                             "rotation_tolerance": 2.0})
placed["events"][0]["placement"]

Constellation figures

project_constellations(wcs, width, height) gives a solved image's constellation stick figures as pixel polylines clipped to it, with a label position for each; constellation_figures() gives every figure and the stars its lines use, for drawing them another way. The figures are the Constellation Lines dataset (Marc van der Sluys, CC BY 4.0): anything that shows them must carry the returned attribution.

Image processing primitives

OpenCV-compatible building blocks from seiza-imgproc, for detection pipelines that need OpenCV's exact numerics without the dependency. All operate on 2D single-channel arrays:

import numpy as np
import seiza

image = np.asarray(..., dtype=np.uint8)      # (height, width)

blurred = seiza.gaussian_blur(image, sigma=1.4)      # uint8 or float32
denoised = seiza.median_blur3(image)
edges = seiza.canny(blurred, low=10, high=80)
binary = seiza.otsu_binary(image)
grown = seiza.dilate(binary, shape="rect", ksize=3)

contours = seiza.find_contours(grown)        # list of (n, 2) int32 arrays
areas = [seiza.contour_area(c) for c in contours]

# Edge-aware smoothing and multi-scale structure removal (float inputs).
flat = seiza.dt_filter(guide, src, sigma_spatial=10.0, sigma_color=30.0)
stars_plus_noise = seiza.remove_structures(image.astype(np.float64), layers=4)

RC-Astro external tools

Drive RC-Astro's standalone rc-astro CLI (BlurXTerminator, StarXTerminator, NoiseXTerminator) on image files. The CLI publishes its own contract — read each tool's parameters from rc_astro_tool_schema rather than hard-coding flags, because they change between CLI builds. Runs need an installed, licensed CLI; the schema's licensed flag and cli_version/ml_version (cache-key material) come from the same probe.

import seiza

print(seiza.rc_astro_locate())               # None when not installed

schema = seiza.rc_astro_tool_schema("sxt")
for p in schema.parameters:
    print(p.name, p.type, p.default, p.min, p.max)

# Star removal, keeping the stars image as a sidecar so starless and
# stars can be stretched independently.
run = seiza.rc_astro_process_file(
    "sxt", "stack.fits", "starless.fits",
    parameters={"stars": True},
    progress=lambda fraction: print(f"{fraction:.0%}"),
)
print(run.primary, run.sidecars, run.device)

A linear stack is far outside the tools' 0..1 range, and they clamp it. rc_astro_process_image runs a tool on an array and returns the result in the input's units. An array within 0..1 keeps its scale; any other is offset by its most negative sample and divided by 65535 or its span, whichever is larger, as PixInsight opens a 16-bit camera frame. NaN and infinite samples, such as a stack's uncovered edges, come back unchanged.

processed = seiza.rc_astro_process_image(tool, stack, parameters=values)
processed.image                            # same shape and units as stack
print(processed.cli_version, processed.ml_version, processed.device)

Ctrl-C cancels a running tool within half a second, as does an optional cancel=lambda: ... predicate returning true. A False switch emits the CLI's --no- negation, overriding a true default. A child completely silent for ten minutes is killed; a first run downloads ML models, which rc-astro download-models handles ahead of time.

PixInsight MLDenoise

Denoise an image with MLDenoise, PixInsight's neural-network noise reduction, through an installed PixInsight with the MachineLearning module. MLDenoise reads its network from a model file (.xmlm) that PixInsight distributes separately. Seiza writes the image to a temporary XISF file, runs PixInsight headless to apply MLDenoise, and reads the result back. On Linux without a display, PixInsight runs under xvfb-run.

import seiza

print(seiza.pixinsight_locate())            # None when not installed

model = "/path/to/MLDenoise_v42.xmlm"
header = seiza.pixinsight_ml_denoise_model(model)
print(header.mono, header.rgb, header.header_sha256)

denoised = seiza.pixinsight_ml_denoise_image(stack, model, amount=0.9)
denoised.image                              # same shape and units as stack
print(denoised.pixinsight_version, denoised.seconds)

seiza.pixinsight_ml_denoise_file(
    "stack.fits", "denoised.fits", model,
    precision="high",
    local_support=seiza.MlDenoiseLocalSupport(clip_low=0.0, background=0.25, smoothness=2),
)

PixInsight works on 0 to 1. An array whose samples lie within 0 to 1 keeps its scale. Any other array is offset by its most negative sample and divided by 65535 or its span, whichever is larger, which is how PixInsight opens a 16-bit camera frame; the result is mapped back. NaN and infinite samples, such as a stack's uncovered edges, come back unchanged. A one-channel image needs a model with a mono route, and an RGB image one with an RGB route: mono and rgb on pixinsight_ml_denoise_model(model) show which routes a model has. The options are PixInsight's, with its ranges and defaults:

MlDenoiseLocalSupport takes clip_low (-20 to 20, in robust sigma from the median), background (0 to 1, the mask's mean background level) and smoothness (0 to 7 pixels, the mask's median-filter radius).

A refused request, such as an option out of range or a model that cannot denoise the image, raises ValueError; a missing model raises FileNotFoundError; a run that fails raises RuntimeError. PixInsight shows license dialogs even in automation mode, and the time limit stops a run that waits on one.

Key a cached result on the model's header_sha256, the pixinsight_version (which includes PixInsight's build) and the options. Two releases of a model can share weights_sha256 and still give different results. PixInsight reports no progress. Ctrl-C cancels a run within a quarter of a second and closes PixInsight, as does an optional cancel=lambda: ... predicate returning true. Runs in one process start PixInsight one at a time, and a run whose PixInsight failed to start because another instance took its slot is started again.

Parallax fly-through videos

Make the video seiza parallax-video makes, from a starless image, its stars and a plate solution: the stars fly at their Gaia distances as the camera moves toward a point of the nebula. rc_astro_split_stars splits a stretched image with StarXTerminator, and solve_blind gives the WCS.

import seiza

starless, stars = seiza.rc_astro_split_stars(image)   # float32 (h, w, 3), 0 to 1
video = seiza.ParallaxVideo(
    starless, stars, solution.wcs,
    focus=(5997, 2351), start="whole", dolly=0.8, seconds=16, size="4k",
    overlay=True, watermark=True,
    progress=lambda event: print(event.get("message", "")),
)
print(video.summary["background_distance_pc"], video.frame_count)

first = video.frame(0)                  # uint8 (height, width, 3)
for frame in video:                     # every frame in order, e.g. for PyAV
    container_stream_encode(frame)
video.render(lambda frame, index: encoder.write(frame), format="bgra")
video.write("iris.mp4")                 # through ffmpeg; encoder="png" for frames

A tour visits several places: tour= takes stops, each a --stop string or a dict. seiza.plan_parallax_tour(width, height, wcs, size=...) plans one from the catalogued objects in the field and returns it for editing:

plan = seiza.plan_parallax_tour(width, height, solution.wcs, size="1080p-portrait")
for stop in plan["tour"]:
    print(stop["name"], stop["focus"], stop["dolly"])
stops = [stop for stop in plan["tour"] if stop["name"] != "LDN 1267"]
video = seiza.ParallaxVideo(starless, stars, solution.wcs, size="1080p-portrait",
                            tour=stops, distance_focus=plan["focus"])

auto_tour=True plans and renders in one go. Each planned stop's title is its target's name; tour_titles=True shows it low in the frame while the camera drifts through the stop, fading in and out, so change or drop titles as you would stops. A tour ends on a longer final drift that slows to rest; tour_loop=True instead ends where it began, so the last frame leads into the first and the video loops. Images with transparency go over the frames too: a stop's card (a path) shows in place of its title, intro_card and outro_card over the start and the end for card_seconds, and watermark={"image": "logo.png", "corner": "top_left"} puts a logo in a corner, with text, size and opacity keys as well.

Depth planes and dust work as in the command and the C API. planes= takes (mask, distance_pc) or (mask, distance_pc, name), each mask a path or a 2-D array, white where the image belongs to the plane; auto_planes=True finds planes from the catalogued objects and the dust; and distance_overrides={"NGC 7822": 900} gives catalogued objects distances of your own. summary["planes"] lists the planes found. dust_map= (a path or a 2-D array, white where the dust is opaque) replaces the dust counted from the stars; video.dust_map() reads the dust back as an array, so it can be painted over, video.with_dust_map(painted) swaps a map in without preparing again, and video.write_dust_map("dust.png") saves it. Once a video is made, video.plan_tour(auto_tour) plans a tour on its scene, each stop framed at the depth of the plane its target lies on, as auto_tour= would; plan_parallax_tour plans from the catalogue alone:

video = seiza.ParallaxVideo(starless, stars, solution.wcs, auto_planes=True)
plan = video.plan_tour({"targets": 6})
toured = video.reconfigure(tour=plan["tour"], tour_titles=True, tour_loop=True)

distance_focus is where the nebula's distance is measured; focus only aims the camera, so moving it never changes the scene's depth. To film the same scene again with another camera, tour, size or quality, call reconfigure. It returns a new video that shares the prepared scene, so nothing is loaded, solved or looked up again, and leaves the first as it was; both can draw frames at once. Options that change the scene, such as the distances, distance_focus, star placement or dust, raise ValueError. summary["fit"] gives each zoom and pan as asked and as used, and each tour stop's dolly, after the camera was fitted to stay inside the image and to show little of the fill behind depth planes, and inside, false when some frames still show past the image's edge. motion_zoom gives the peak extra magnification fitted during a transition when both endpoint views already fit; the endpoint lenses stay unchanged. Held spins can still increase their own stop's fitted zoom. For stops with pull_back_to_whole_image=True, pull-back centering uses overall zoom multiplied by the stop's zoom:

wide = video.reconfigure(size="4k", quality="high", rotate_deg=(0, -360))
print(wide.summary["fit"]["zoom"])      # (asked, used)

An image of another shape than the frame is cut to fill it by default. framing="fit" shows all of it at the widest, with black beside it, and a size ending -image takes the image's own shape: "4k-image" films a 4:3 image at 3840x2880, all of it. crop=(x, y, width, height) films part of the image, and every other image point, such as focus and tour stops, counts in the cropped image; pass the same crop to plan_parallax_tour:

video = seiza.ParallaxVideo(starless, stars, solution.wcs,
                            crop=(800, 400, 6000, 4000), size="4k-image")
square = video.reconfigure(size=(1080, 1080), framing="fit", start="whole")

Images may be paths (PNG, JPEG, TIFF, FITS or XISF) or arrays, float32 from 0 to 1 or uint8; 16-bit or float files make the best videos, as 8 bits flatten bright star cores and band smooth nebula. Options are the command's, in snake case, with frame sizes such as "720p", "1080p", "4k", "1080p-portrait" or "4k-image". Distances come from the star distance file when installed, or the archives unless online=False. render stops when its callback returns False; Ctrl-C or cancel stop preparing, reconfigure, render and write.

Predicted satellite tracks

After a solve, predict which satellites crossed the image while the shutter was open. Predictions come from orbital elements — they are never pixel detections. The exposure must be one continuous shutter-open interval (not a stack's total integration) and needs an observer location:

sats = seiza.SatelliteCatalog.fetch_celestrak()   # cached; ~2h refresh floor
# or offline / historical: seiza.SatelliteCatalog.open("elements.json")

result = sats.tracks_in_footprint(
    solution.wcs, width, height,
    start="2026-07-19T06:12:00Z",     # Unix seconds, RFC 3339, or tz-aware datetime
    duration_s=120.0,                 # or end=, in the same forms as start
    latitude=42.466, longitude=-71.1516, altitude_m=150.0,
)
for track in result.tracks:           # highest elevation first
    print(track.label, track.max_elevation_deg, track.clipped_segments)

Element records older than seven days are reported in result.stale_elements and skipped rather than silently extrapolated (max_element_age_s=None overrides). CelesTrak rate-limits repeated downloads: keep reusing one cache directory, and check sats.cache_state and sats.warning after fetch_celestrak().

Working on the bindings

This directory is its own cargo workspace, so cargo fmt --all and cargo clippy --workspace at the repository root never reach it. Run its checks from here, which is what CI does:

cd seiza-py
cargo fmt --check
cargo clippy --all-targets -- -D warnings
maturin develop && python -m pytest tests/ -q

Notes

• Solving and detection release the GIL; other Python threads keep running.
• Parallel work runs in Rayon's global pool, which by default has a thread for every core. To cap it, set RAYON_NUM_THREADS before the first call into seiza, for example in the environment that starts Python.
• Catalog files are memory-mapped and SHA-256 verified at download time; fetch_catalogs caches under the platform cache directory (override with cache_dir= or SEIZA_CACHE_DIR).
• seiza.StarCatalog.from_stars([...]) builds a small in-memory catalog for tests and synthetic fields.

License

Apache-2.0

Subtract stellar continuum before color composition

fit_continuum(stars, sources, options=None) fits premeasured, isolated common stars. measure_and_fit_continuum(narrowband, references, sources, positions=None, fit_options=None, photometry_options=None) measures them on aligned, star-containing linear mono float32 arrays. Omit positions to detect stars. Each function returns the core JSON: either {"status": "fitted", "plan": ...} or {"status": "skipped", "reason": ..., "detail": ..., "diagnostics": ...}. A skipped fit has no plan; decide how to recover before composing an image.

source = lambda name, revision: dict(name=name, revision=revision,
                                    transfer="linear", saturation=65535)
sources = dict(width=ha.shape[1], height=ha.shape[0],
               narrowband=source("ha", "ha:registered:v1"),
               references=[source("red", "red:registered:v1")])
fit = seiza.measure_and_fit_continuum(ha, [red], sources,
                                    fit_options={"minimum_stars": 12})
if fit["status"] == "fitted":
    emission = seiza.subtract_continuum(ha, [red], sources, fit["plan"])

saturation is the physical saturation level in each plane's working units, not its observed maximum. Use 1 for samples scaled from a sensor saturation level to one. Change revision whenever pixels, scale, background treatment, registration or masks change. Names must be distinct; the one to three references and their order define coefficient order. Plans reject changed sources, scales (through revisions), grids, order or non-linear transfer. Python option dicts accept snake_case; returned plan/diagnostic keys match the C ABI's camelCase JSON. Input arrays are copied before releasing the GIL.

The default fit requires 12 common unsaturated stars, SNR 20, width ratios at most 1.3, centroid differences at most one pixel, and training/held-out RMS fractional errors at most 0.15. Every fourth accepted star is held out; at least 70% of both sets must remain after rejection. Aperture/sky radii default to 8/12/20 pixels. fit_options accepts all ContinuumFitOptions fields; photometry_options accepts aperture_radius, annulus_inner, annulus_outer. Manual plans have schemaVersion=1, the same sources, a nonnegative coefficients list, origin="manual", and diagnostics=None. subtract_continuum(narrowband, references, sources, plan) returns a new signed array N - sum(a_i * B_i), keeping negative values and NaN for missing contributing coverage. It never clips, rescales, stretches or corrects filter line transmission. A broadband reference can contain emission itself, so the result is an emission residual rather than calibrated pure line flux.

Use producer clipping masks and add emission

When no physical working-unit saturation ceiling is known, set each source's clipping_revision and call measure_and_fit_continuum_masked(narrowband, references, sources, clipping_masks, *, positions=None, fit_options=None, photometry_options=None). Supply a C-contiguous uint8 mono mask for each input, narrowband first: 0 is usable and 1 is clipped. The masks select stars; they do not alter pixels or signed subtraction. Automatic measurement examines up to 4096 candidates independently of the usable-star budget, applies eligibility checks before that budget, samples across stellar brightness, and checks the fit against the complete eligible pool. eligibleStars, sampledStars and samplingStabilityRms report these decisions. A failed broader validation returns a skipped outcome. Manual subtraction can omit saturation evidence.

blend_emission_rgb(broadband, residuals, lines, *, negative="preserve") projects signed mono residuals into linear RGB. Each line contains its exact plan, reference_channels (RGB indices in reference order), an artistic amount from 0 to 5, and optional explicit background_offset (default 0). The original residual remains signed. negative="positive_contribution" limits only the added contribution. Fitted plans record referenceFluxRatios; manual plans with multiple references must supply these measured or explicit ratios. For stellar ratios r_j=F_Bj/F_B0, projection uses k_j=sum(a_ir_i)/r_j and adds amount(E-background_offset)/k_j. Independent input gains commute with this conversion. It is an empirical emission residual, not pure line flux.

Calibrate the broadband stars before adding emission. To apply their measured gains to an already normalized RGB contribution, call add_rgb_contribution(base, contribution, *, gains=(1,1,1)); offsets apply to the base only. blend_display_structure(narrowband, broadband, amount) mixes prepared starless RGB display images in [0,1], before restoring RGB stars once. Amount 0 keeps narrowband, 1 keeps broadband, and intermediate amounts require common coverage. All three functions return new arrays and preserve inputs.

Validate profiles and export calibrated line flux

measure_and_fit_continuum and its masked twin accept photometry_options={"psf_matching": {"target_fwhm": 5.0}} for explicit Gaussian broadening on copies. Omit target_fwhm to choose 1.05 times the broadest measured profile. The plan records raw/prepared revisions and kernels; subtract_continuum replays them on the original arrays. Saturation evidence is still required for every original input. Held-out pixel validation defaults to max_pixel_residual=0.15 and remains mandatory for automatic host use. The result's pixelValidation report distinguishes cores/halos from aperture flux error. Gaussian matching can still fail this check.

balance_narrowband(..., star_mask_radius=8) excludes star circles from curve statistics while preserving image coverage. None includes stars. reduce_green_excess(image, amount) returns a copy of bounded display RGB, reducing green above both red and blue. Use it after palette mapping and before restoring an RGB-star layer. The input is unchanged.

solve_emission_flux(residuals, plans, calibration) returns (images, report). The calibration mapping supplies revision, unit, lines, exact sources, narrowband_response, reference_response and optional residual_covariance. Rows follow the plan/line order, references follow each plan. Responses include all exposure/gain/throughput factors in the exact input units. The operation subtracts the fitted broadband responses before solving coupled lines. Signed samples and common coverage are retained. The report embeds calibration/plans and propagated covariance. Missing covariance means unknown uncertainty; response, coefficient and zero-level uncertainty are not included. Filter labels cannot supply a physical flux calibration.

Matching also requires target_fwhm <= aperture_radius and annulus_inner >= 1.5 * target_fwhm. The default 8/12/20 photometry therefore permits targets up to 8 pixels. Larger targets require larger photometry radii in the options JSON (with a hard target limit of 32 pixels and outer sky-annulus radius at most 128 pixels). This keeps the target profile inside the measurement aperture and away from its sky annulus.

Calibrated export also requires explicit residualZeroLevels (Python alias residual_zero_levels), one additive zero level per residual in its exact sample units. Zero is an explicit calibration declaration. The solver uses M F = E − residualZeroLevels; it never guesses these values from a sky median or the display baseline. Background preparation may retain an additive reference level, so response calibration alone cannot establish a physical zero point. Zero-level uncertainty is outside the propagated residual covariance unless the producer includes it there.

Image-based validation also audits isolated candidates excluded for clipping or saturation. They never train the fit. Their observed residual profiles must all be measurable and below the same absolute limit; clippedStars, clippedProfiles and clippedMaxAbsoluteFraction record this separate check. It prevents a clean unsaturated sample from hiding clipped-core artifacts and does not infer lost flux. The held-out report keeps both its 90th-percentile and maximum profile error.

screen_rgb_stars_reduced(nebula, original, starless, reduction) returns display RGB with a reduced screen-star layer. Reduction 0 preserves ordinary screen recombination; 1 uses full one-pixel cross erosion of star intensity. A shared RGB gain preserves star-layer color ratios, while invalid neighbours/borders retain the unreduced layer. The inputs and starless nebula are unchanged. This changes display prominence/apparent size, not linear or PSF measurements.

Calibrated export requires continuumCoefficients (Python alias continuum_coefficients): exact per-plan coefficient snapshots in reference order. Residual zero-level calibration depends on these weights. A changed fit/manual weight is rejected even when the raw source revisions are unchanged; update the explicit zero levels and snapshots after recalibration. Reports and cache identities retain the snapshots.

stretch_rgb_stars(original, starless, amount=5) takes an aligned linear RGB original and its linear starless result, extracts the residual, and returns stars stretched with the SetiAstro-style rational tone. Amount [0,8] controls prominence: lower is fainter; 0 keeps the normalized residual without nonlinear lift. Shared source normalization handles physical units. This is the tone formula, without PixInsight's optional saturation/SCNR processes. screen_rgb_star_layer(nebula, stars, reduction=0) screens the display layer over a display nebula with optional linked star reduction [0,1]. Inputs are unchanged; missing coverage masks whole RGB pixels. Run star separation before stretching.

Joint or reference canvas

channels = [seiza.CanvasChannel(ha_lights, flat=ha_flat),
            seiza.CanvasChannel(oiii_lights, flat=oiii_flat)]
ha = seiza.LiveStacker.open_fits_joint_canvas(channels, channel=0, options=options)
for path in ha_lights:
    ha.push_fits(path)
ha.save_context("union.seiza-stack")
oiii = ha.open_fits_on_canvas(oiii_lights, flat=oiii_flat, options=options)
for path in oiii_lights:
    oiii.push_fits(path)

Joint mode builds the union of connected filter footprints; the first group's first light anchors geometry. Each filter keeps its own calibration and normalization, and missing coverage stays NaN. CanvasChannel accepts paths and optional bias, dark, flat, dark_exposure_seconds, exposed as read-only properties. open_fits_joint_canvas also accepts max_canvas_megasamples and cancel like open_fits_expanding; channel is zero-based (default 0).

For faster reference mode when a saved canvas already covers all filters, use LiveStacker.open_context(...).open_fits_on_canvas(...). Joint mode adds a geometry pass over every group; do it once and reuse its result. Stellar overlap is required between panels; normalization needs overlap within each filter.

Expanded mosaics with calibration sessions

LiveStacker.open_fits_expanding_sessions prepares every original light with its assigned calibration before planning the canvas. It returns an empty stack; push every original file, including the anchor, in its original order:

sessions = [
    {"bias": "bias-a.fits", "flat": "flat-a.fits"},
    {"bias": "bias-b.fits", "flat": "flat-b.fits",
     "tolerances": {"light_rotation": {"basis": "measured", "degrees": 90.0}}},
]
stacker = seiza.LiveStacker.open_fits_expanding_sessions(
    paths, sessions, [0, 1, 0], max_canvas_megasamples=256)
for path in paths:
    stacker.push_fits(path)

Each session supplies a complete gate; omitted values use native defaults. A fitted-pedestal recipe supplies flat (already normalized) and pedestal, without bias or dark. normalized_full_scale=65535 rescales only images that declare normalized 0:1 bounds, matching PSF Guard's loader. Assignments are immutable and survive checkpoints and reintegration. Changed sources or masters are refused. See session planning for file identity, memory and pipeline behavior.

LiveStacker.open_fits_sessions_with_reference_image(reference_image, paths, sessions, assignments, channel_reference) takes the same sessions on one reference grid instead of a canvas. paths[channel_reference], calibrated by its own session, sets normalization and noise weighting, and the stack starts empty. Stacks that share the reference image, the channel reference and its session put their lights on one grid, scale and weight, so each session can be stacked on its own and the stacks combined, as distributed stacking does.

Distributed stacking

seiza.run_channel(plan, store) stacks a channel as partial stacks, one per session, into a store, and combines them, running only the steps whose outputs the store lacks (see distributed stacking). Each session aligns to the reference image and is normalized and weighted against the reference light, calibrated by its own recipe, so the partial stacks combine into what one stack of every light gives, apart from rejection, which runs within each session.

reference = {"image": "n1/L_001.fits", "light": "n1/L_001.fits",
             "calibration": {"bias": "n1/bias.fits"}}
sessions = [{"lights": n1_lights, "calibration": {"bias": "n1/bias.fits"}},
            {"lights": n2_lights, "calibration": {"bias": "n2/bias.fits"}}]
outcome = seiza.run_channel(
    {"reference": reference, "sessions": sessions,
     "options": {"reintegration": {}, "contributor": "me"}},
    "store")
tile = seiza.PartialStack.open(outcome["tile"]["path"])
tile.mean, tile.weight, tile.coverage   # numpy arrays

seiza.SessionStep.plan(reference, session, options) computes one session's key from its files' contents; run() stacks it, and run_into(store) returns the stored partial stack or stacks and stores it. With "keep_context": True in the options, a session that gains lights is extended from its stored context, and the answer's extended_from names the step it was extended from; LiveStacker.add_session_lights(paths, session) is the call underneath for a host that pushes lights itself. PartialStack.read_identity reads a file's identity cards without its pixels, and PartialStack.combine combines stacks on one grid.

Calibration library: grading and selection

Python holds the calibration library as Rust hosts such as PSF Guard do (seiza_calibration::library): the same objects and calls, no request round trip. Records are snake_case dicts of the Calibration* types (frames, lights, quality rows); nothing here opens a file.

rules = seiza.CalibrationRules(
    {"dark_choice": "complete"},                              # the server's layer
    {rig_uuid: {"name": "SpaceCat61", "values": {"dark_reach_days": 400}}},
)
nights = seiza.CalibrationNights(night_boundaries, frames=frames, rules=rules)
quality = seiza.CalibrationQuality(quality_rows)
grades = quality.regrade(frames, rules, nights)   # current measurements only
frames = quality.keep_out(frames, grades)         # each frame's kept-out reason
library = seiza.CalibrationLibrary(frames, groups=seiza.FlatGroups(grouping_rows), rules=rules)
for light in lights:
    selection = library.select(light, rules.for_light(light))
    session = seiza.calibration_resolve_session(selection, quality)
    print(session.key, session.recipes)

CalibrationRules(global_settings, rigs) checks every value the way a settings panel saves it and raises ValueError naming the setting (and rig). for_rig and for_light give a CalibrationRigRules: value(key), source(key), set_by(key), mode(check), notes(), tolerances(), grouping_limits(), and admits_sensor(light, frame), the sensor test selection asks. With allow_flat_gain_offset_mismatch on, it lets a flat master whose dark_subtracted is set differ from the light in gain and offset. calibration_settings_schema() lists the settings in the order a panel shows them, with labels, help, ranges and defaults, and calibration_setting_check(key, value) judges one value: valid, the error in words, the value as text, and whether it means the default.

CalibrationNights holds where each frame's night begins; given the frames and rules, each rig's night_zone setting applies over its frames. CalibrationQuality.regrade grades only current, fingerprint-matching measurements and returns CalibrationGrades; measurement_state(frame) tells stale, failed and missing measurements apart. keep_out reads each frame's kept_out as the stray-light mark an older build left, so give it the frames as the catalog stores them, not frames it returned. A FlatGroups book is built once from stored grouping rows, with each rig's limits from rules.for_rig(rig).grouping_limits(), and reused across selections.

select(light, rules, flat_association="automatic", dark_association="automatic") returns a CalibrationSelection: each kind's frames nearest first, the kept-out matches, the flat basis and family, and match_tolerances(). Set an association to "one_set" to take every eligible frame of that kind as one supplied set, without the automatic 64-input cap; check the set first with library.check_associations(rules, ...), which raises ValueError for frames that cannot integrate together. CalibrationSelection(frames, rules, ...) builds the selection of frames a host chose or verified itself.

calibration_resolve_session(selection, quality, flat_star_masking=False, flat_stability=None) gives the CalibrationSession: its key (policy version 5: what its masters are made from, as a set, so lights whose masters take the same frames share one), the stack gate tolerances, and the bias, dark-flat, dark and flat recipes. An unverified selection is a preview: verify the sources and resolve again before building. Lights with equal keys calibrate identically; the first one stands for the session.

Shared calibration master execution

Build a session's recipes in order, telling the session what happened: session.record_failure(kind, reason) blocks every master built on a failed one. Once every master has built or failed, session.plan_flat_use(lights, without_masters="fit_pedestal") decides how the flat applies to the session's own lights; with neither a bias nor a dark master it asks for a pedestal fitted from application["sample_lights"], which session.accept_pedestal(estimates, pinned=None) accepts or refuses (without_masters="force" applies the flat regardless). calibration_coherent_cohort(kind, frames, rules) gives the frames one master integrates and those it drops, for previews. A host that builds its own masters, as PSF Guard does, makes the decisions a session's recipes are made of with the same calls: calibration_adopted_master (the external master a light takes), calibration_dark_level_screening, calibration_pedestal_sample_lights, calibration_pedestal_hint, calibration_pedestal_from_fits, calibration_family_key, calibration_pedestal_moved, calibration_flat_set_label and calibration_describe_finding.

FlatRetry(first, candidates, rules, session=...) decides a flat master's retries from measured build feedback (FlatAttempt): leave a drifting minority out, else take the next flat session, else keep the first set with a note, or withhold it when the policy keeps no unstable set. The session gives the retry policy and held-back flats (a supplied set never takes another session), and rules come from the selection; without a session, pass held_back and policy:

retry = seiza.FlatRetry(first_attempt, raw_flats, selection.rules, session=session)
decision = retry.next()
while decision["action"] in ("try_trimmed", "try_alternate"):
    built = build(decision["inputs"])              # None when no master came of it
    (retry.trimmed if decision["action"] == "try_trimmed" else retry.alternate)(built)
    decision = retry.next()

rebuild_first says the first build is a memo to build again for real. A complete host execution loop shows retries, memo rebuilding, publication and cleanup.

master_flat_stability(master=..., inputs=[{path, exposure_seconds, skipped}], bias=None, dark=None, dark_exposure_seconds=None, normalized_full_scale=None) measures each source against a master file. It returns algorithm_version (2), feature_pixels, drift_threshold, disagreement and indexed errors. Missing or incompatible inputs count as drifting; skipped inputs have zero disagreement. normalized_full_scale=65535 preserves PSF Guard's declared-unit XISF handling. Mono, RGB and CFA data use separate sample planes. See the execution contract.

master_flat_stability_images(master, inputs, *, bias=None, dark=None, input_metadata=None, bayer=None, dark_exposure_seconds=None, dark_bias_subtracted=False) measures cached float32 NumPy images without file I/O. Arrays must be C-contiguous mono (H,W) or interleaved RGB (H,W,3); raw CFA uses mono plus bayer={pattern, x_offset, y_offset}. All pixels already use the same units. Metadata is omitted/empty or aligned to inputs, with exposure_seconds, bayer and skipped. The bindings copy buffers for the call; the host retains its original caches. File and cached-image measurement execute the same algorithm and return FlatSetStability.

Analysis baselines and stellar SNR

The result dictionary uses snake_case (median_star_snr, median_normalized_noise). Frames now include median unsaturated-star peak SNR. Filter groups expose measured-frame medians of star count, HFR, FWHM, eccentricity, noise, background and stellar SNR, plus gain-normalized noise. Raw noise/background retain native input units; normalized noise uses reference units. Missing measurements are null/None. See measurement definitions for population, units and limitations. Thresholds and rejection remain caller-owned.

Completed-stack reports

Capture context at the same point as the result, before finishing or adding more frames. Reporting reads the result's existing maps; it does not reread images or run star detection.

context = stacker.report_context()
snapshot = stacker.reintegrate()  # or stacker.snapshot()
report = snapshot.report(context)
Path("stack.html").write_text(seiza.stack_report_html(report), encoding="utf-8")
Path("stack.json").write_text(json.dumps(report, indent=2), encoding="utf-8")

DrizzleResult.report(context) uses accumulated drop area times frame weight, not frame counts. Final sample decisions refer to the reintegration grid. Both report methods accept no context when input history is unavailable. Use context from the same run; matching counts alone cannot establish identity. Restored contexts can lack exposures and registration measurements. In-memory inputs without file provenance have unknown exposures. The admitted ledger does not retain failed attempts; preserve push outcomes separately if needed. The standalone HTML has light/dark themes, bounded per-channel support maps, expandable input records, and a download of its embedded JSON.

Combining completed stacks

Report generation is optional. A collection contains independent stack results; it never re-registers images or infers shared sky coverage. Labels are nonempty and unique (ignoring case and surrounding whitespace), and input order determines the initial selection. Missing filters and different output sizes are supported. RGB components remain inside their own stack, distinct from collection labels.

The collection JSON is { "schema_version": 1, "title": "SHO", "stacks": [{ "label": "S", "report": <PostStackReport> }, ...] }. The HTML embeds the validated collection for download, uses no network resources, and supports keyboard controls, search, 100-row pagination and printing all stacks. Rejected attempts appear when supplied in host provenance. Empty context remains unknown, not zero exposure. Maps are bounded summaries; this is not a crop tool.

Rendering saved JSON does not reread images or compute checksums. Building the original post-stack report scans existing support maps and copies ledger metadata; CLI provenance additionally hashes input/output files. Reporting never enables reintegration by itself.

Python: html = seiza.stack_collection_html(collection) returns a string and raises ValueError for invalid collections or nested reports.

Cached flat evidence and correction workflow

The measuring stages of flat grouping are the ones PSF Guard runs, on records the host keeps; the host reads headers and pixels when asked.

• calibration_flat_sessions(frames, nights, rules) gives the raw flat sessions to measure, one night of one family each. A FlatSessionPlan names the combined() flats to read, its response_key(), and the pending_row(pedestal) or measured_row(pedestal, read, bin, response, unread) to keep, response being flat_build_response's answer. PedestalTable(frames, quality) gives each session's pedestal (for_session) and the pedestal under lights (light_pedestal).
• calibration_family_questions(library, rules, rig, lights_of, header) asks the lights of each imaging session which response family their flats come from. lights_of(session) returns a session's lights as {path, at, angle, filter, target}; header(path) returns a light's header record, or None. A FamilyQuestion that is answered needs nothing read. Otherwise read its sample() lights binned, take each representatives() session's response plane with its pedestal out (flat_response_without_pedestal), and keep decide(representatives, binned, planes, pedestal)'s row in the grouping's family_evidence.
• calibration_change_question(groups, event, session, lights) asks the lights shot during a change which side of it they stand on, for each (event, session) of groups.sessions_in_windows(); those in groups.sessions_needing_evidence() have no answer yet. For each ChangeFilter, measure change_map(before, after) from its map_responses() and keep ChangeQuestion.row(fits) with each filter's fit(binned, reference, map, pedestal) in the grouping's evidence.

Store the rows unchanged and build the FlatGroups book again. No lights read means retry later; lights read but no useful background or fit are kept as undecided. The core keeps PSF Guard's nearest-angle tie policy and refuses missing fits as proof of similarity.

Correction masters use flat_change_map, preserving the full compact map; flat_evidence_change_map masks border vignetting specifically for light fitting. Apply a correction with flat_apply_change_map(..., bin=..., divide=...). Keep the original master metadata and include the correction key in cache identity. File verification, atomic publication and failed cache recovery remain host responsibilities.

Manual group_overrides and light_choices are already core features on grouping rows. Feed PSF Guard's existing rows directly; the library selector enforces hard physical gates and resolves stale/orphaned choices. Manual light-choice entries can use an exact full path in file_name. It takes precedence over basename entries for that light, keeping same-named files from different folders separate.

No report, cheap report, rich report

Skip context/report construction entirely when reporting is disabled. The shared prepare_stack_report(report=None, *, level="cheap", inputs=()) helper returns None for level="none" and does not open inputs. Cheap uses the existing completed-stack report; rich adds an explicit source-analysis pass:

report = snapshot.report(context)
cheap = seiza.prepare_stack_report(report, level="cheap")
rich = seiza.prepare_stack_report(report, level="rich", inputs=light_paths)
Path("rich.html").write_text(seiza.stack_report_html(rich), encoding="utf-8")

Rich requires explicit source paths and retains failed measurements. It reads original source pixels without stack calibration, measures stellar quality, and computes sampled equal-weight depth curves; these are not final-stack noise/SNR estimates. Results and per-stage timings are in report["provenance"]["reporting"]["analysis"]. No mode changes the stack. Cheap chart rendering also works on previously saved reports. Downgrading rich to cheap removes the rich-analysis payload without reading source images. The C twin is seiza_prepare_stack_report_json, taking {report, level, inputs} with camelCase report keys and returning camelCase JSON (null for none). Whole-file checksums are CLI provenance behavior, not performed by this helper.

OSC / RGB analysis (schema 3)

Analysis automatically recognizes mono, raw Bayer and three-plane RGB sources. No extra option is needed, including for rich post-stack reports. The input list is read once per source; measurement and registration run per component. max_frames limits source paths; samples is the budget per component.

A schema-3 frames record represents one source/component, not necessarily one file. CFA produces R, G_red (green on a red row), G_blue (green on a blue row), and B; RGB produces R, G, B; mono produces Mono. Source-level read failures or filter exclusions produce one record without sampling metadata. Groups separate (filter, image_kind, component), so CFA/RGB/mono never share normalization or noise baselines. Do not sum exposure across color components.

Frame sampling records image_kind, component, source_width, source_height, spacing_pixels, offset_pixels, declared bayer_pattern, storage_bayer_pattern (row-order adjusted, before offsets), bayer_offsets, and declared row_order. C JSON uses camelCase spellings of these same keys. Native CFA sample (x,y) is at original stored-image pixel (2*x + offset_pixels[0], 2*y + offset_pixels[1]). Each lattice retains every native sample, including odd edges. No demosaicing, green averaging, white balance or new calibration occurs. Unrecognized Bayer declarations and invalid origin offsets are reported as failures rather than guessed as mono.

HFR, FWHM and registration RMS are in original-image pixels (CFA component values multiplied by two); this does not recover details lost to sparse CFA sampling. Undersampled stars can have biased shapes. Noise/background retain native component units; depth noise uses that component's reference units. Registration and photometric normalization are independent per component. A component with too few stars/common patches can keep measurements while its noise curve is unavailable. RGB retains existing demosaicing/resampling correlations; it must not be interpreted as independent raw detector samples.

Both HTML renderers select filter/layout/component explicitly. Schema 1/2 reports remain renderable; new reports use schema/algorithm version 3. Stacking and stretching numerical output versions do not change.

File-backed calibration coordinator

CalibrationFiles(lights, cache, *, library=None, provided=None, rules=None, incomplete="stop", flat_association="automatic", dark_association="automatic", flat_star_masking=False, flat_stability=None, choices=None, group_overrides=None, light_choices=None) uses the same native file coordinator as the desktop C ABI. rules is a CalibrationRules object. analyze(cancel=None, progress=None) grades/groups/matches without building masters; prepare(review_key, cancel=None, progress=None) revalidates the review and publishes each unique dependency and corrected flat atomically. Reports are snake_case records including snapshot, flat_work, missing, light_paths, staged paths, file_sessions and sources. CalibrationFiles.verify(sources) refuses changed inputs/masters before use. Callback exceptions and Ctrl-C are preserved; work releases the GIL.

The coordinator is optional. Database-backed applications can keep using the individual native calibration, quality, session and evidence classes with their own caches and pixels. See file calibration for native Rust/C examples, cache identities, cancellation and source rules.