From 0aa37cc503dccf866159c88985a323fc99537b28 Mon Sep 17 00:00:00 2001 From: laurence Date: Tue, 21 Jul 2026 19:39:30 +0100 Subject: [PATCH] Add the measure and register stages, generic across sessions measure.py produces the numbers used to choose a registration reference and to weight the stack - sky, noise, seeing - and caches each frame's star list, because registration and the plate solve both need it and detection costs far more than reading a small array back. Frames are opened one at a time; a 4096x4096 float32 frame is 67 MB and a session holds ninety of them. register.py aligns everything onto a single reference and combines per filter. One reference for ALL filters, not one per filter, which is what makes the masters pixel-aligned so the colour composite needs no further registration. The reference is the sharpest frame of the filter with the most signal, because the reference sets the output grid and a poor choice costs resolution everywhere, permanently. The case worth the care is NGC 6744: luminance at bin1 4096x4096 and colour at bin2 2048x2048, so frames share neither shape nor pixel scale. Asterism matching already returns a similarity transform including scale, so the maths was never the problem; the trap is assuming the warp output is the same shape as its input, which would write a bin2 frame onto a bin2 grid that silently fails to line up with a bin1 reference. Every warp is now given the reference shape explicitly. Verified rather than assumed: after stacking, the upsampled bin2 colour masters align to the bin1 luminance master to within 0.04 pixels, across about 245 matched stars per channel. Sigma clipping is skipped below three frames, where there is nothing to reject against and clipping would only discard signal - which matters because one session has a single frame per filter. Tested on three sessions covering the awkward shapes: RGB without luminance, SHO with one frame per filter, and the mixed bin1/bin2 LRGB. --- pipeline/measure.py | 127 ++++++++++++++++++++++++++++++++++ pipeline/register.py | 159 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 286 insertions(+) create mode 100644 pipeline/measure.py create mode 100644 pipeline/register.py diff --git a/pipeline/measure.py b/pipeline/measure.py new file mode 100644 index 0000000..2593ee1 --- /dev/null +++ b/pipeline/measure.py @@ -0,0 +1,127 @@ +"""Stage 2: measure every frame and cache its star list. + +Two jobs. It produces the quality numbers used to choose a registration +reference and to weight the stack - sky level, noise, star sharpness - and it +caches each frame's detected stars, because both registration and the plate +solve need them and detection costs far more than reading a small array back. + +Frames are opened one at a time. A 4096x4096 float32 frame is 67 MB and a +session can hold ninety of them; there is no reason to have more than one in +memory at once. +""" +import os + +import numpy as np +import sep + +from astropy.io import fits + +import layout + +MAX_STARS = 500 # cached per frame, brightest first +DETECT_SIGMA = 5.0 + + +def measure_frame(path, max_stars=MAX_STARS): + """Sky, noise, seeing and a star list for one frame.""" + with fits.open(path, memmap=False) as hd: + data = hd[0].data.astype(np.float32) + + bkg = sep.Background(data, bw=64, bh=64, fw=3, fh=3) + sky = float(np.median(bkg.back())) + rms = float(bkg.globalrms) + sub = data - bkg.back() + + objs = sep.extract(sub, DETECT_SIGMA, err=rms, minarea=9, + deblend_cont=0.005) + # Saturated cores and cosmic-ray hits both make poor registration anchors + # and poor seeing estimates, so they are excluded before anything is + # measured from them. + objs = objs[(objs["flag"] == 0) & (objs["npix"] > 12) & + (objs["npix"] < 3000)] + if len(objs) == 0: + del data, sub, bkg + return dict(sky=sky, rms=rms, fwhm=np.nan, nstars=0, + xy=np.zeros((0, 2)), flux=np.zeros(0)) + + order = np.argsort(objs["flux"])[::-1][:max_stars] + objs = objs[order] + radius, _ = sep.flux_radius(sub, objs["x"], objs["y"], 6.0 * objs["a"], + 0.5, normflux=objs["flux"]) + fwhm = float(np.median(radius) * 2.0) + + result = dict( + sky=sky, rms=rms, fwhm=fwhm, nstars=int(len(objs)), + xy=np.column_stack([objs["x"], objs["y"]]).astype(np.float64), + flux=objs["flux"].astype(np.float64), + ) + del data, sub, bkg, objs + return result + + +def run(session, force=False, verbose=True): + """Measure every frame in a session, caching to intermediates/.""" + cache_path = os.path.join(session.root, layout.INTERMEDIATES, + "_measure.npz") + os.makedirs(os.path.dirname(cache_path), exist_ok=True) + + cached = {} + if os.path.exists(cache_path) and not force: + z = np.load(cache_path, allow_pickle=True) + cached = {k: z[k] for k in z.files} + + store, rows = {}, [] + for frame in session.frames: + key = frame.name + if f"{key}|xy" in cached and not force: + xy = cached[f"{key}|xy"] + stats = cached[f"{key}|stats"].item() + else: + stats = measure_frame(frame.path) + xy = stats.pop("xy") + flux = stats.pop("flux") + store[f"{key}|flux"] = flux + store[f"{key}|xy"] = xy + store[f"{key}|stats"] = np.array(stats, dtype=object) + rows.append((frame, stats)) + if verbose: + print(f" {frame.filter:10s} {frame.name[-28:]:28s} " + f"sky={stats['sky']:8.1f} rms={stats['rms']:7.2f} " + f"fwhm={stats['fwhm']:5.2f}px stars={stats['nstars']:4d}") + + # Carry forward anything already cached that was not re-measured. + for k, v in cached.items(): + store.setdefault(k, v) + np.savez_compressed(cache_path, **store) + return rows + + +def load(session): + """Star lists and stats for a session, as measured earlier.""" + path = os.path.join(session.root, layout.INTERMEDIATES, "_measure.npz") + if not os.path.exists(path): + raise FileNotFoundError( + f"{path} not found - run the measure stage first") + z = np.load(path, allow_pickle=True) + xy = {k.split("|")[0]: z[k] for k in z.files if k.endswith("|xy")} + stats = {k.split("|")[0]: z[k].item() for k in z.files + if k.endswith("|stats")} + return xy, stats + + +def choose_reference(session, stats): + """Pick the frame everything else is aligned to. + + The reference sets the output grid, so it should be the sharpest frame of + the filter with the most signal - which is also usually the filter with the + finest sampling. Choosing a poor reference costs resolution in every other + frame, permanently, because registration can only resample onto it. + """ + best_filter = session.filters[0] + candidates = [f for f in session.by_filter(best_filter) + if stats.get(f.name, {}).get("nstars", 0) >= 8] + if not candidates: + candidates = session.by_filter(best_filter) + ref = min(candidates, + key=lambda f: stats.get(f.name, {}).get("fwhm", np.inf)) + return ref diff --git a/pipeline/register.py b/pipeline/register.py new file mode 100644 index 0000000..c4a522d --- /dev/null +++ b/pipeline/register.py @@ -0,0 +1,159 @@ +"""Stage 3: register every frame onto one grid and combine per filter. + +Registration matches asterisms between the cached star lists rather than +cross-correlating pixels: matching a few hundred coordinates is far cheaper +than comparing 15 megapixel frames, and it copes with the field rotating +between the east and west sides of the meridian. + +**One reference for all filters, not one per filter.** That is what makes the +masters pixel-aligned, so the colour composite needs no further registration. + +The awkward case this has to survive: NGC 6744 was shot with luminance at bin1 +(4096x4096) and colour at bin2 (2048x2048). Frames do not share a shape or a +pixel scale. Because astroalign returns a similarity transform - rotation, +translation AND scale - the maths already handles it; what does not handle it +is assuming the output is the same shape as the input. Every warp is therefore +given the reference's shape explicitly, and the bin2 frames are resampled up +onto the bin1 grid. +""" +import os + +import astroalign as aa +import numpy as np +from astropy.io import fits +from astropy.stats import sigma_clip +from skimage.transform import warp + +import layout +import measure as measure_mod + +aa.MIN_MATCHES_FRACTION = 0.6 +aa.NUM_NEAREST_NEIGHBORS = 8 + + +def combine(cube, weights, sigma=3.0): + """Weighted sigma-clipped mean, in row blocks to bound peak memory. + + sigma_clip allocates a mask and float64 intermediates; on a full cube that + can triple peak usage, which matters on a machine with a few GB free. + """ + n, ny, nx = cube.shape + out = np.zeros((ny, nx), np.float32) + w = np.asarray(weights, np.float32)[:, None, None] + step = 256 + for y0 in range(0, ny, step): + y1 = min(y0 + step, ny) + block = cube[:, y0:y1, :] + if n >= 3: + clipped = sigma_clip(block, sigma=sigma, maxiters=1, axis=0, + masked=True, copy=True) + good = ~clipped.mask + else: + # With one or two frames there is nothing to reject against, and + # clipping would just throw away signal. + good = np.ones(block.shape, bool) + finite = np.isfinite(block) + good &= finite + vals = np.where(finite, block, 0.0) + wb = np.broadcast_to(w, block.shape) * good + denom = wb.sum(axis=0) + denom[denom == 0] = np.nan + out[y0:y1, :] = np.nansum(vals * wb, axis=0) / denom + del block, finite, good, vals, wb, denom + return np.nan_to_num(out, nan=0.0) + + +def run(session, verbose=True): + """Register and stack every filter; write masters to stacks/masters/.""" + xy, stats = measure_mod.load(session) + ref = measure_mod.choose_reference(session, stats) + ref_xy = xy[ref.name] + ref_shape = ref.shape + if verbose: + print(f" reference {ref.name[-32:]} ({ref.filter}, " + f"fwhm {stats[ref.name]['fwhm']:.2f} px, " + f"{ref_shape[1]}x{ref_shape[0]})") + + out_dir = os.path.join(session.root, layout.MASTERS) + os.makedirs(out_dir, exist_ok=True) + written = {} + + for filt in session.filters: + frames = session.by_filter(filt) + planes, weights, used = [], [], [] + header0 = None + for frame in frames: + with fits.open(frame.path, memmap=False) as hd: + data = hd[0].data.astype(np.float32) + if header0 is None: + header0 = hd[0].header.copy() + if frame.name == ref.name: + reg = data + else: + src = xy.get(frame.name) + if src is None or len(src) < 8: + print(f" {frame.name[-28:]}: too few stars - " + f"dropped") + del data + continue + try: + tform, _ = aa.find_transform(src, ref_xy) + except Exception as exc: # noqa: BLE001 + print(f" {frame.name[-28:]}: registration failed " + f"({type(exc).__name__}) - dropped") + del data + continue + # output_shape is the whole point: without it a bin2 frame + # would be written onto a bin2-sized grid and silently fail to + # line up with a bin1 reference. + reg = warp(data, inverse_map=tform.inverse, + output_shape=ref_shape, order=3, mode="constant", + cval=np.nan, preserve_range=True).astype(np.float32) + del data + sky = float(np.nanmedian(reg)) + reg -= sky + planes.append(reg) + rms = stats.get(frame.name, {}).get("rms", 1.0) or 1.0 + weights.append(1.0 / (rms * rms)) + used.append(frame) + + if not planes: + print(f" {filt}: no frames survived registration - skipped") + continue + + cube = np.stack(planes, axis=0) + del planes + # Normalise for transparency: scale each frame so its bright signal + # matches the group, so a frame through thin cloud cannot drag the + # mean down. + levels = np.array([np.nanpercentile(p, 99.5) for p in cube]) + ref_level = np.nanmedian(levels) + for i, lv in enumerate(levels): + if lv > 0: + cube[i] *= float(ref_level / lv) + + master = combine(cube, weights) + nframes = cube.shape[0] + total_exp = sum(f.exptime for f in used) + del cube + + hdr = header0 + for key in ("CBLACK", "CWHITE", "PEDESTAL", "HISTORY"): + hdr.remove(key, ignore_missing=True, remove_all=True) + hdr["FILTER"] = filt + hdr["NCOMBINE"] = (nframes, "frames in this master") + hdr["EXPTOTAL"] = (total_exp, "[s] total integration") + hdr["STACKREF"] = (ref.name[:60], "registration reference frame") + hdr["STACKALG"] = ("sigma-clipped weighted mean" if nframes >= 3 + else "weighted mean (too few frames to clip)", + "combine method") + hdr["IMAGETYP"] = "Master Light" + path = os.path.join(out_dir, f"master-{filt}.fit") + fits.PrimaryHDU(master.astype(np.float32), hdr).writeto( + path, overwrite=True) + written[filt] = path + if verbose: + print(f" {filt:10s} {nframes:3d} frames " + f"{total_exp/60:6.1f} min -> master-{filt}.fit") + del master + return written, ref