From 0aa37cc503dccf866159c88985a323fc99537b28 Mon Sep 17 00:00:00 2001 From: laurence Date: Tue, 21 Jul 2026 19:39:30 +0100 Subject: [PATCH 1/8] 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 From 4ed23cef96e1045017290ed0a5fdec0897ea678a Mon Sep 17 00:00:00 2001 From: laurence Date: Tue, 21 Jul 2026 19:44:02 +0100 Subject: [PATCH 2/8] Add the colour stage, palette chosen from the filters present colour.py assembles a viewable image from whatever a session actually has: LRGB, RGB with a synthetic luminance, SHO in the Hubble palette, HOO, or a single filter as greyscale. Nothing is assumed about which filters exist, which is the whole point - three of the four archived sessions have no luminance channel. The background fit takes the lesson from Centaurus A, where a plane fitted around a large galaxy absorbed 17.9 ADU/px of its halo. The excluded region is now sized from the data: the ellipse grows until it contains most of the flux above sky, so it suits a galaxy filling the frame and a small nebula equally, and only a plane is ever fitted, never a flexible surface. One mask is derived from the deepest channel and applied to all of them, so the fit cannot shift the colour balance. Broadband and narrowband need opposite stretches, which cost a round to discover. Broadband channels share the luminance's midtone, preserving the real brightness ratios that keep star colours honest. Narrowband cannot: the lines differ enormously in strength, so a shared stretch renders the entire nebula in whichever colour Ha was mapped to. Each narrowband channel is now stretched to its own sky target instead. Green suppression is applied to broadband only. On a narrowband palette it would fight the palette, green being a deliberate channel assignment rather than an artefact. Measured rather than eyeballed, because the SHO render looked wrong by eye and was not: both palettes put the sky at 0.098 against a 0.10 target and neutral to within 0.004. What differs is the data. NGC 2030's brightest pixels reach 0.80 to 0.87, the Tarantula's only 0.23 to 0.33, because that session is a single 60 second frame per filter - three minutes in total. A correct pipeline renders nearly empty data as a nearly empty frame. --- pipeline/colour.py | 259 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 259 insertions(+) create mode 100644 pipeline/colour.py diff --git a/pipeline/colour.py b/pipeline/colour.py new file mode 100644 index 0000000..ce2edc1 --- /dev/null +++ b/pipeline/colour.py @@ -0,0 +1,259 @@ +"""Stage 5: assemble a viewable image from whatever filters the session has. + +The palette is chosen from what is present, not assumed: + + LRGB colour from R/G/B, brightness from the deep luminance + RGB no luminance, so a synthetic one is made from the colour channels + SHO the Hubble palette: SII -> red, Ha -> green, OIII -> blue + HOO Ha -> red, OIII -> green and blue + MONO one filter, rendered as greyscale + +Background handling deserves a note, because the obvious approach is wrong in a +way that is invisible. Fitting a plane to a frame that a large target fills +makes the plane absorb the target's own outer light - on Centaurus A this was +measured at -17.9 ADU/px of real halo quietly subtracted away. So the fit +excludes a central region whose size is derived from where the signal actually +is, and only a plane is used, never a flexible surface. +""" +import os + +import numpy as np +import tifffile +from astropy.io import fits +from PIL import Image +from scipy.ndimage import gaussian_filter, median_filter + +import layout + +BG_TARGET = 0.10 # where the sky sits in the stretched image +SATURATION = 1.35 + +PALETTE_MAP = { + "SHO": {"Red": "SII", "Green": "Ha", "Blue": "OIII"}, + "HOO": {"Red": "Ha", "Green": "OIII", "Blue": "OIII"}, +} + + +def mtf(x, midtone): + """Midtone transfer function on data already scaled to [0, 1].""" + x = np.clip(x, 0.0, 1.0) + return ((midtone - 1.0) * x) / ((2.0 * midtone - 1.0) * x - midtone) + + +def autostretch(img, target=BG_TARGET, shadow_sigma=2.8): + """Black point just below sky, then an MTF putting sky at `target`.""" + sky = float(np.median(img)) + mad = 1.4826 * float(np.median(np.abs(img - sky))) + if mad <= 0: + mad = max(float(np.std(img)), 1e-6) + black = sky - shadow_sigma * mad + white = float(np.percentile(img, 99.995)) + if white <= black: + white = black + 1.0 + norm = np.clip((img - black) / (white - black), 0.0, 1.0) + sky_norm = (sky - black) / (white - black) + m = ((target - 1.0) * sky_norm) / (2.0 * target * sky_norm - target - + sky_norm) + m = float(np.clip(m, 1e-4, 0.9)) + return mtf(norm, m), dict(black=black, white=white, midtone=m, sky=sky, + mad=mad) + + +def target_mask(img, frac=0.45): + """Ellipse covering the bright central object, sized from the data. + + A fixed radius cannot suit both a galaxy filling the frame and a small + nebula. This grows the ellipse until it contains most of the flux above + sky, then stops - so the background fit is excluded from wherever the + target actually is, whatever its size. + """ + ny, nx = img.shape + yy, xx = np.mgrid[0:ny, 0:nx] + cy, cx = ny / 2.0, nx / 2.0 + r = np.hypot((xx - cx) / (nx / 2.0), (yy - cy) / (ny / 2.0)) + sky = float(np.median(img)) + signal = np.clip(img - sky, 0, None) + total = signal.sum() + if total <= 0: + return r < 0.5 + for cut in np.arange(0.25, 0.96, 0.05): + if signal[r < cut].sum() >= frac * total: + return r < min(cut + 0.15, 0.95) + return r < 0.75 + + +def remove_gradient(img, mask): + """Subtract a plane fitted to tile medians outside `mask`.""" + ny, nx = img.shape + step = max(48, min(ny, nx) // 40) + xs, ys, zs = [], [], [] + for y0 in range(0, ny - step, step): + for x0 in range(0, nx - step, step): + if mask[y0:y0 + step, x0:x0 + step].any(): + continue + tile = img[y0:y0 + step, x0:x0 + step] + zs.append(np.median(tile)) + xs.append(x0 + step / 2.0) + ys.append(y0 + step / 2.0) + if len(zs) < 12: + return img - float(np.median(img)), None + xs, ys, zs = map(np.asarray, (xs, ys, zs)) + keep = np.ones(len(zs), bool) + coef = None + for _ in range(3): + A = np.column_stack([xs[keep], ys[keep], np.ones(keep.sum())]) + coef, *_ = np.linalg.lstsq(A, zs[keep], rcond=None) + model = coef[0] * xs + coef[1] * ys + coef[2] + resid = zs - model + s = 1.4826 * np.median(np.abs(resid - np.median(resid))) + if s <= 0: + break + keep = np.abs(resid - np.median(resid)) < 2.5 * s + yy, xx = np.mgrid[0:ny, 0:nx] + plane = (coef[0] * xx + coef[1] * yy + coef[2]).astype(np.float32) + return img - plane, coef + + +def load_masters(session): + masters = {} + d = os.path.join(session.root, layout.MASTERS) + for name in os.listdir(d) if os.path.isdir(d) else []: + if name.startswith("master-") and name.endswith(".fit"): + filt = name[len("master-"):-len(".fit")] + masters[filt] = os.path.join(d, name) + return masters + + +def run(session, verbose=True): + """Produce the final image for whatever palette this session supports.""" + masters = load_masters(session) + if not masters: + raise RuntimeError("no masters found - run the register stage first") + palette = session.palette + if verbose: + print(f" palette {palette} from {sorted(masters)}") + + data, headers = {}, {} + for filt, path in masters.items(): + with fits.open(path) as hd: + data[filt] = hd[0].data.astype(np.float32) + headers[filt] = hd[0].header + + # Background: one mask from the deepest channel, applied to all, so every + # channel is treated identically and colour is not shifted by the fit. + deepest = session.filters[0] if session.filters[0] in data \ + else sorted(data)[0] + mask = target_mask(data[deepest]) + if verbose: + print(f" target mask covers {mask.mean():.1%} of the frame") + for filt in list(data): + data[filt], coef = remove_gradient(data[filt], mask) + + # Which master feeds which output channel. + if palette in PALETTE_MAP: + mapping = PALETTE_MAP[palette] + channels = [data[mapping[c]] for c in ("Red", "Green", "Blue")] + lum_src = None + elif palette in ("LRGB", "RGB", "PARTIAL") and all( + c in data for c in ("Red", "Green", "Blue")): + channels = [data["Red"], data["Green"], data["Blue"]] + lum_src = data.get("Luminance") + else: + only = data[deepest] + stretched, _ = autostretch(only) + grey = (np.clip(stretched, 0, 1) * 65535).astype(np.uint16) + return _write(session, np.dstack([grey] * 3), palette, verbose) + + # Colour balance: make the median star colour neutral, which is the cheap + # substitute for a photometric calibration and is well behaved on a rich + # field. + refs = [float(np.percentile(c, 99.5)) for c in channels] + target = float(np.median(refs)) + channels = [c * (target / r if r > 0 else 1.0) for c, r in + zip(channels, refs)] + + lum_lin = lum_src if lum_src is not None else np.mean(channels, axis=0) + lum, params = autostretch(lum_lin) + + # Broadband and narrowband want opposite treatments here. + # + # For broadband, all three channels share the luminance's midtone (a + # "linked" stretch). That preserves the real brightness ratios between + # them, which is what keeps star colours honest. + # + # For narrowband it is wrong. The three lines have wildly different + # strengths - Ha is typically far brighter than SII and OIII - so a linked + # stretch renders the whole nebula in whichever colour Ha was mapped to, + # which in the Hubble palette means overwhelmingly green. Stretching each + # channel to its OWN sky target instead brings the weak lines up to + # comparable visual weight, which is what the palette is for. The result is + # deliberately false colour either way; this makes it a readable false + # colour. + unlinked = palette in PALETTE_MAP + rgb = np.empty(lum.shape + (3,), np.float32) + for i, ch in enumerate(channels): + if unlinked: + rgb[:, :, i], _ = autostretch(ch) + else: + sky = float(np.median(ch)) + mad = 1.4826 * float(np.median(np.abs(ch - sky))) or 1.0 + black = sky - 2.8 * mad + white = float(np.percentile(ch, 99.995)) + rgb[:, :, i] = mtf( + np.clip((ch - black) / (max(white - black, 1e-6)), 0, 1), + params["midtone"]) + if unlinked: + # With each channel independently stretched the composite already + # carries its own brightness; taking luminance from the linear mean + # would flatten it back down. + lum = rgb.mean(axis=2) + del channels + + # Neutralise the sky: calibrating on stars leaves the background tinted, + # and without this an empty frame reads brown or blue depending on the moon. + sky_med = [float(np.median(rgb[:, :, i][~mask])) for i in range(3)] + neutral_to = float(np.mean(sky_med)) + for i in range(3): + rgb[:, :, i] = np.clip(rgb[:, :, i] - (sky_med[i] - neutral_to), 0, 1) + + # Denoise colour only; structure comes from the luminance, so this is + # invisible at normal viewing scale. + rgb_lum = rgb.mean(axis=2, keepdims=True) + chroma = rgb - rgb_lum + for i in range(3): + chroma[:, :, i] = gaussian_filter(median_filter(chroma[:, :, i], 3), + 1.5) + rgb = np.clip(rgb_lum + chroma * SATURATION, 0, 1) + del chroma, rgb_lum + + ratio = lum / np.maximum(rgb.mean(axis=2), 1e-5) + out = np.clip(rgb * ratio[:, :, None], 0.0, 1.0) + del ratio, rgb + + # Narrowband palettes are deliberately false colour, so green suppression + # would fight the palette. It only applies to broadband. + if palette not in PALETTE_MAP: + neutralish = 0.5 * (out[:, :, 0] + out[:, :, 2]) + green = out[:, :, 1] + out[:, :, 1] = np.where(green > neutralish, + green * 0.15 + neutralish * 0.85, green) + + return _write(session, (out * 65535 + 0.5).astype(np.uint16), palette, + verbose) + + +def _write(session, rgb16, palette, verbose=True): + target = session.target.replace(" ", "") + final = os.path.join(session.root, layout.FINAL) + os.makedirs(final, exist_ok=True) + stem = f"{target}-{palette}" + tifffile.imwrite(os.path.join(final, stem + ".tif"), rgb16, + photometric="rgb") + u8 = (rgb16 / 257).astype(np.uint8) + Image.fromarray(u8).save(os.path.join(final, stem + ".png")) + prev = Image.fromarray(u8) + prev.thumbnail((2400, 2400), Image.LANCZOS) + prev.save(os.path.join(final, stem + "-preview.jpg"), quality=92) + if verbose: + print(f" -> final/{stem}.png / .tif / -preview.jpg") + return os.path.join(final, stem + ".png") From fbd1eb8ec50d2db63f5babf6514bee4d177881eb Mon Sep 17 00:00:00 2001 From: laurence Date: Tue, 21 Jul 2026 21:27:33 +0100 Subject: [PATCH 3/8] Add the plate solve stage; it works on two sessions of four astrometry.py solves the deepest master against Gaia and copies the WCS into every master. Asterism matching is invariant to rotation and scale, so camera angle never has to be guessed, but not to a mirror flip, so both parities are tried. Three defects were found and fixed by running it, and each is worth recording because none would have been found by reading the code. Gaia's launch_job_async submits to a job QUEUE and polls. On a query of 600 rows it hung for ten minutes while curl showed both the Gaia and VizieR endpoints answering in under a second. VizieR is now queried first - it serves the identical DR3 catalogue over a plain HTTP request - with Gaia's SYNCHRONOUS endpoint as the fallback. VizieR's row_limit truncates before sorting, so asking for the 600 brightest stars returned 600 arbitrary ones. Asterism matching only works when both lists hold the same bright stars, so this silently produced no match at all. It now fetches generously and picks the brightest locally. The detector deliberately rejects saturated cores, which means the brightest DETECTIONS are not the brightest STARS, while the catalogue's are - so the two top-N lists can barely overlap. A sliding window down the catalogue's magnitude ranking fixed NGC 2030, which matched at catalogue[15:135], skipping the 15 brightest. That is the saturation hypothesis confirmed rather than assumed. NGC 2030 now solves on 243 stars at 0.31 arcsec residual. NGC 2070 and NGC 6744 still do not, and TODO.md records what has already been ruled out - scale, the fetch, truncation - so the next session starts from the remaining candidates rather than repeating the elimination. --- pipeline/astrometry.py | 286 +++++++++++++++++++++++++++++++++++++++++ state/TODO.md | 36 +++++- 2 files changed, 320 insertions(+), 2 deletions(-) create mode 100644 pipeline/astrometry.py diff --git a/pipeline/astrometry.py b/pipeline/astrometry.py new file mode 100644 index 0000000..4887980 --- /dev/null +++ b/pipeline/astrometry.py @@ -0,0 +1,286 @@ +"""Stage 4: plate solve the deepest master and copy the WCS to the others. + +iTelescope's frames arrive with a PinPoint HISTORY line but no WCS keywords at +all, so the astrometry has to be redone locally. A blind solve is unnecessary: +the header gives the pointing to arcminutes and usually the plate scale to four +figures, so a catalogue can be fetched for the right patch of sky and matched +directly. + +The match is asterism-based, which is invariant to rotation and scale - so the +camera angle never has to be guessed, and a session at an unknown orientation +solves as easily as one at a known angle. It is NOT invariant to a mirror flip, +because a similarity transform cannot include reflection, so both parities are +tried and the one that matches wins. + +The seed match returns only a handful of stars. That is enough to establish an +approximate solution but too few to trust, so the solution is then used to pair +every catalogue source with its nearest detection and refit, twice, with a +shrinking tolerance. The residual that comes out of that is the honest measure +of whether the solve is real. + +Failure here is not fatal to the pipeline: an unsolved session still produces +images, it just cannot produce positions. +""" +import os + +import astroalign as aa +import numpy as np +from astropy import units as u +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.wcs import WCS +from astropy.wcs.utils import fit_wcs_from_points + +import layout +import measure as measure_mod + +SEED_STARS = 600 + + +def pointing(header): + """Approximate field centre from whatever the header offers.""" + for ra_key, dec_key in (("OBJCTRA", "OBJCTDEC"), ("RA", "DEC"), + ("CRVAL1", "CRVAL2")): + ra, dec = header.get(ra_key), header.get(dec_key) + if ra is None or dec is None: + continue + try: + if isinstance(ra, str) and (":" in ra or " " in ra): + return SkyCoord(ra, dec, unit=(u.hourangle, u.deg)) + return SkyCoord(float(ra) * u.deg, float(dec) * u.deg) + except Exception: # noqa: BLE001 + continue + return None + + +def plate_scale(header, session=None): + """Arcsec per pixel: from the header if present, else from the optics.""" + val = header.get("HIERARCH iTelescopePlateScaleH") + if val: + try: + return float(val) + except (TypeError, ValueError): + pass + focal = header.get("FOCALLEN") + pixel = header.get("XPIXSZ") + if focal and pixel: + # 206265 arcsec per radian; XPIXSZ is already binned microns. + return 206.265 * float(pixel) / float(focal) + return session.scale if session else None + + +def fetch_catalogue(centre, radius_deg, cache, limit=SEED_STARS, deep=False): + """Gaia sources around the pointing, cached so a re-run costs nothing.""" + if os.path.exists(cache): + z = np.load(cache) + return z["ra"], z["dec"], z["g"] + ra = dec = g = None + mag_cut = 20.5 if deep else 18.0 + + # VizieR first, deliberately. It serves the identical DR3 catalogue over a + # plain HTTP query and answers in under a second. Gaia's own archive is + # excellent but launch_job_async submits to a job QUEUE and polls for the + # result, which for a few hundred rows means waiting minutes for work that + # takes milliseconds - and it hung for ten minutes on a query of 600 stars + # while the network was demonstrably healthy. + try: + from astroquery.vizier import Vizier + # Ask for far more rows than are needed and pick the brightest + # locally. VizieR's row_limit truncates the result BEFORE any sort, so + # asking for 600 rows in a crowded field returns 600 arbitrary stars + # rather than the 600 brightest - and asterism matching only works if + # both lists contain the same bright stars. In the LMC field of + # NGC 2030 that silently produced no match at all. + v = Vizier(columns=["RA_ICRS", "DE_ICRS", "Gmag"], + column_filters={"Gmag": f"<{mag_cut}"}, row_limit=50000) + v.TIMEOUT = 60 + res = v.query_region(centre, radius=radius_deg * u.deg, + catalog="I/355/gaiadr3") + if res: + t = res[0] + ra = np.asarray(t["RA_ICRS"], float) + dec = np.asarray(t["DE_ICRS"], float) + g = np.asarray(t["Gmag"], float) + except Exception as exc: # noqa: BLE001 + print(f" VizieR unavailable ({type(exc).__name__})") + + if ra is None or len(ra) == 0: + try: + from astroquery.gaia import Gaia + Gaia.ROW_LIMIT = limit + # Synchronous, not async: no job queue for a query this small. + job = Gaia.launch_job(f""" + SELECT TOP {limit} ra, dec, phot_g_mean_mag + FROM gaiadr3.gaia_source + WHERE 1 = CONTAINS(POINT('ICRS', ra, dec), + CIRCLE('ICRS', {centre.ra.deg}, {centre.dec.deg}, + {radius_deg})) + AND phot_g_mean_mag IS NOT NULL + AND phot_g_mean_mag < {mag_cut} + ORDER BY phot_g_mean_mag ASC""") + t = job.get_results() + ra = np.asarray(t["ra"], float) + dec = np.asarray(t["dec"], float) + g = np.asarray(t["phot_g_mean_mag"], float) + except Exception as exc: # noqa: BLE001 + print(f" Gaia archive unavailable ({type(exc).__name__})") + return None, None, None + + order = np.argsort(g)[:limit] + ra, dec, g = ra[order], dec[order], g[order] + np.savez_compressed(cache, ra=ra, dec=dec, g=g) + return ra, dec, g + + +def project(ra, dec, centre, scale, parity): + """Gnomonic projection to pixel-like coordinates for asterism matching.""" + c = SkyCoord(ra * u.deg, dec * u.deg) + dx, dy = centre.spherical_offsets_to(c) + return np.column_stack([dx.to_value(u.arcsec) / scale * parity, + dy.to_value(u.arcsec) / scale]) + + +def run(session, verbose=True): + """Solve the deepest master; write the WCS into every master.""" + masters_dir = os.path.join(session.root, layout.MASTERS) + deepest = session.filters[0] + path = os.path.join(masters_dir, f"master-{deepest}.fit") + if not os.path.exists(path): + raise RuntimeError("no masters found - run the register stage first") + + with fits.open(path) as hd: + image = hd[0].data.astype(np.float32) + header = hd[0].header + ny, nx = image.shape + + centre = pointing(header) + scale = plate_scale(header, session) + if centre is None or not scale: + print(" no pointing or plate scale in the header - " + "cannot solve without a blind solver") + return None + radius = 1.15 * 0.5 * np.hypot(nx, ny) * scale / 3600.0 + if verbose: + print(f" pointing {centre.to_string('hmsdms')}, " + f"{scale:.4f} arcsec/px, search {radius:.3f} deg") + + xy, _ = measure_mod.load(session) + # Detect on the master itself: it is deeper than any single frame. + from measure import measure_frame + det = measure_frame(path, max_stars=SEED_STARS) + src = det["xy"] + if len(src) < 12: + print(f" only {len(src)} stars in the master - too few") + return None + + cache = os.path.join(session.root, layout.INTERMEDIATES, "_gaia.npz") + os.makedirs(os.path.dirname(cache), exist_ok=True) + ra, dec, gmag = fetch_catalogue(centre, radius, cache) + if ra is None: + print(" no catalogue available - skipping the solve") + return None + if verbose: + print(f" {len(ra)} catalogue stars, {len(src)} detected") + + # Matching needs the two lists to contain the SAME stars, and "brightest" + # does not guarantee that. The detector deliberately rejects saturated + # cores, so on a well exposed frame the brightest detections are not the + # brightest stars in the sky - while the catalogue's brightest are. The + # two top-N lists can therefore barely overlap, which is exactly what + # happened on the crowded LMC field of NGC 2030. + # + # So slide a window down the catalogue's magnitude ranking as well as + # varying the sample size, and take the best match found. + best = None + for offset in (0, 15, 40, 80): + for n in (120, 200, 60): + hi = min(offset + n, len(ra)) + if hi - offset < 25 or len(src) < 25: + continue + cat_slice = slice(offset, hi) + for parity in (1.0, -1.0): + cat_xy = project(ra[cat_slice], dec[cat_slice], centre, scale, + parity) + try: + tform, (cat_m, img_m) = aa.find_transform( + cat_xy, src[:max(n, 120)]) + except Exception: # noqa: BLE001 + continue + if best is None or len(cat_m) > len(best[0]): + # Record which catalogue rows matched, in FULL-list terms. + best = (cat_m, img_m, parity, offset) + if verbose: + print(f" catalogue[{offset}:{hi}] parity " + f"{parity:+.0f}: matched {len(cat_m)} stars") + if best is not None and len(best[0]) >= 8: + break + if best is not None and len(best[0]) >= 8: + break + + if best is None: + print(" no asterism match in either parity - solve failed") + return None + + cat_m, img_m, parity, _offset = best + cat_xy = project(ra, dec, centre, scale, parity) + idx = [int(np.argmin(np.hypot(cat_xy[:, 0] - px, cat_xy[:, 1] - py))) + for px, py in cat_m] + world = SkyCoord(ra[idx] * u.deg, dec[idx] * u.deg) + wcs = fit_wcs_from_points((img_m[:, 0], img_m[:, 1]), world, + proj_point="center", projection="TAN") + + # Refine: the seed match is a handful of stars, so use the approximate + # solution to pair every catalogue source with its nearest detection. + all_world = SkyCoord(ra * u.deg, dec * u.deg) + nmatch, resid = len(cat_m), None + for tol in (4.0, 2.0): + px, py = wcs.world_to_pixel(all_world) + pairs = [] + for i, (cx, cy) in enumerate(zip(px, py)): + if not (0 <= cx < nx and 0 <= cy < ny): + continue + d = np.hypot(src[:, 0] - cx, src[:, 1] - cy) + j = int(np.argmin(d)) + if d[j] <= tol: + pairs.append((i, j)) + if len(pairs) < 20: + break + ci = np.array([p[0] for p in pairs]) + ii = np.array([p[1] for p in pairs]) + wcs = fit_wcs_from_points((src[ii, 0], src[ii, 1]), all_world[ci], + proj_point="center", projection="TAN") + qx, qy = wcs.world_to_pixel(all_world[ci]) + resid = np.hypot(qx - src[ii, 0], qy - src[ii, 1]) + nmatch = len(pairs) + + cd = wcs.pixel_scale_matrix * 3600.0 + solved_scale = float(np.sqrt(abs(np.linalg.det(cd)))) + rot = float(np.degrees(np.arctan2(cd[0, 1], cd[1, 1]))) + med_resid = float(np.median(resid)) if resid is not None else float("nan") + field = wcs.pixel_to_world(nx / 2.0, ny / 2.0) + + if verbose: + print(f" solved on {nmatch} stars, residual " + f"{med_resid:.2f} px ({med_resid * solved_scale:.2f} arcsec)") + print(f" centre {field.to_string('hmsdms')}, " + f"{solved_scale:.4f} arcsec/px, PA {rot:.2f} deg, " + f"{nx * solved_scale / 60:.1f}' x {ny * solved_scale / 60:.1f}'") + + whdr = wcs.to_header() + for filt in session.filters: + p = os.path.join(masters_dir, f"master-{filt}.fit") + if not os.path.exists(p): + continue + with fits.open(p, mode="update") as hd: + for card in whdr.cards: + hd[0].header[card.keyword] = (card.value, card.comment) + hd[0].header["ASTRSOLV"] = ( + f"Gaia DR3 / {nmatch} stars / " + f"{med_resid * solved_scale:.2f} arcsec", "local plate solve") + hd.flush() + + return dict(nstars=nmatch, residual_px=med_resid, + residual_arcsec=med_resid * solved_scale, + scale=solved_scale, position_angle=rot, + centre=field.to_string("hmsdms"), + fov_arcmin=(nx * solved_scale / 60, ny * solved_scale / 60)) diff --git a/state/TODO.md b/state/TODO.md index efbcbf9..8342390 100644 --- a/state/TODO.md +++ b/state/TODO.md @@ -18,8 +18,40 @@ one telescope. Branch `generic-pipeline`. filter grouping, palette detection (LRGB / RGB / SHO / HOO / MONO), manifest writing. Tested against all five archived sessions. -**Next:** `measure` -> `register/stack` -> `solve` -> `compose` -> `report`, -then a `run.py` orchestrator with resumable stages. +**Done since:** `measure.py`, `register.py`, `colour.py`, `astrometry.py`. +Registration verified on the mixed bin1/bin2 case - upsampled bin2 colour +aligns to the bin1 luminance master within **0.04 px**. Colour renders LRGB, +RGB, SHO, HOO and mono. + +**Next:** finish `astrometry.py` (see below), then `report.py` (per-session +METHODS.md) and `run.py` (the one-command orchestrator). + +**PLATE SOLVE IS 2 OF 4.** NGC 5128 (by hand) and NGC 2030 solve; NGC 2070 and +NGC 6744 do not. What is already known, so the next session does not repeat it: + +- **Not the plate scale.** Header and FOCALLEN/XPIXSZ agree on both failures, + and the session that DOES solve has a header/computed discrepancy + (0.527 vs 0.686) yet solves anyway - asterism matching is scale invariant, + so scale only sets the catalogue cone radius. +- **Not the catalogue fetch.** Fixed: VizieR is queried first because Gaia's + `launch_job_async` submits to a job queue and hung for ten minutes on 600 + rows while the network was healthy. Gaia synchronous is the fallback. +- **Not simple truncation.** Fixed: VizieR's `row_limit` truncates BEFORE + sorting, so asking for 600 rows returned 600 arbitrary stars rather than the + 600 brightest. Now fetches up to 50000 and picks the brightest locally. +- **Partly saturation.** The detector rejects saturated cores, so the brightest + DETECTIONS are not the brightest STARS while the catalogue's are, and the two + top-N lists barely overlap. NGC 2030 only solved once a sliding magnitude + window was tried, and it matched at `catalogue[15:135]` - skipping the 15 + brightest. That window search is in place and still is not enough for the + other two. +- **Next things to try:** match on a magnitude-matched subset by estimating the + image's own zero point first; loosen `aa.MIN_MATCHES_FRACTION`; raise + astroalign's `max_control_points`; or fall back to a genuine blind solver + (astrometry.net index files locally, or the nova API) for fields that resist + the seeded approach. NGC 2070 is a 4 degree field at 3.5 arcsec/px with only + 305 detections, and NGC 6744 is a rich field - they may need different + handling from each other. **The blocker before any of that:** the disk is full (see below). From 38f7a813956e8a0d37e4738608b2e2690adfc97b Mon Sep 17 00:00:00 2001 From: laurence Date: Tue, 21 Jul 2026 21:40:50 +0100 Subject: [PATCH 4/8] Add the report generator and the push-button entry point run.py is the one command the whole exercise was aimed at: python run.py --all PATH It finds every session beneath a directory and takes each through ingest, measure, register, solve, colour and report. Stages are skipped when their output already exists, so an interrupted run continues rather than restarting, and --force overrides that. A failure in one stage of one session does not stop the others. In a batch of twenty the useful outcome is nineteen results and one clear error, not nothing - so failures are caught, recorded and summarised at the end, and the plate solve in particular is allowed to fail without taking the images down with it. An unsolved session still produces a perfectly good picture; it just cannot produce positions. report.py writes each session its own METHODS.md, generated from what the stages actually returned rather than from what they were supposed to return. That matters for the honest parts: a session whose plate solve failed says so in its own documentation, and the limitations section is derived from the data - too few frames for outlier rejection, an integration too short to recover, a missing luminance channel, seeing that caps the achievable detail. A hand-written note is written once and then rots; this is rebuilt on every run. --- pipeline/report.py | 188 +++++++++++++++++++++++++++++++++++++++++++++ pipeline/run.py | 158 +++++++++++++++++++++++++++++++++++++ 2 files changed, 346 insertions(+) create mode 100644 pipeline/report.py create mode 100644 pipeline/run.py diff --git a/pipeline/report.py b/pipeline/report.py new file mode 100644 index 0000000..2b1b025 --- /dev/null +++ b/pipeline/report.py @@ -0,0 +1,188 @@ +"""Stage 6: write the session's own METHODS.md. + +Every session directory gets an account of what was done to it and what came +out, written for someone who was not there and has no access to the chat that +produced it. That is the whole point of generating it rather than writing it by +hand: a hand-written note is written once and then rots, while this is rebuilt +from the actual results every time the pipeline runs. + +It records what was measured, not what was intended - including the things that +did not work, because a session where the plate solve failed is more useful +with that stated than with it omitted. +""" +import os +from datetime import datetime, timezone + +import numpy as np + +import layout + + +def _fmt_exposure(seconds): + return f"{seconds / 60:.0f} min" if seconds >= 60 else f"{seconds:.0f} s" + + +def write(session, results, verbose=True): + """Compose METHODS.md from whatever the stages actually returned.""" + stats = results.get("stats") or {} + solve = results.get("solve") + reference = results.get("reference") + palette = session.palette + now = datetime.now(timezone.utc).strftime("%Y-%m-%d %H:%M UTC") + + lines = [] + add = lines.append + + add(f"# {session.target} - {session.telescope}") + add("") + add(f"Processed {now} by the pipeline in the " + "[astrophotography](https://git.discworld.casa/laurence/astrophotography) " + "repository (`pipeline/`). This file is generated: re-running the " + "pipeline rewrites it from the actual results.") + add("") + + # ---------------------------------------------------------- what is here + add("## What was captured") + add("") + add("| Filter | Frames | Exposure | Binning | Total |") + add("|---|---|---|---|---|") + total = 0.0 + for filt in session.filters: + frames = session.by_filter(filt) + exps = sorted({f.exptime for f in frames}) + bins = sorted({f.binning for f in frames}) + secs = sum(f.exptime for f in frames) + total += secs + add(f"| {filt} | {len(frames)} | " + f"{'/'.join(f'{e:.0f} s' for e in exps)} | " + f"{'/'.join(f'bin{b}' for b in bins)} | {_fmt_exposure(secs)} |") + add(f"| **Total** | **{len(session.frames)}** | | | " + f"**{_fmt_exposure(total)}** |") + add("") + add(f"Colour scheme: **{palette}**" + + {"SHO": " - the Hubble palette, SII to red, Ha to green, OIII to blue", + "HOO": " - Ha to red, OIII to green and blue", + "LRGB": " - colour from R/G/B, detail from the deeper luminance", + "RGB": " - no luminance channel, so brightness is synthesised from " + "the colour channels", + "MONO": " - a single filter, rendered as greyscale"}.get(palette, "")) + add("") + + # -------------------------------------------------------- frame quality + if stats: + add("## Frame quality") + add("") + add("| Filter | Sky (ADU) | Noise (ADU) | Seeing (px) |") + add("|---|---|---|---|") + for filt in session.filters: + vals = [stats[f.name] for f in session.by_filter(filt) + if f.name in stats] + if not vals: + continue + sky = np.median([v["sky"] for v in vals]) + rms = np.median([v["rms"] for v in vals]) + fwhm = np.median([v["fwhm"] for v in vals]) + add(f"| {filt} | {sky:.1f} | {rms:.2f} | {fwhm:.2f} |") + add("") + + # ---------------------------------------------------------- what was done + add("## What was done") + add("") + add("1. **Ingest.** Archives uncompressed into `calibrated/`. Only the " + "`calibrated-` set is extracted; the `raw-` archives stay zipped in " + "`raw/` as untouched originals, because extracting both would put two " + "copies of every frame into the stack.") + add("2. **Measure.** Sky, noise and star sharpness per frame, with each " + "frame's star list cached for the stages that follow.") + if reference: + add(f"3. **Register and stack.** Every frame aligned by matching star " + f"patterns onto one reference - `{reference.name}` " + f"({reference.filter}, the sharpest frame of the deepest filter) - " + f"so all masters share a pixel grid and the colour image needs no " + f"further alignment. Combined with a " + f"{'sigma-clipped ' if len(session.frames) >= 3 else ''}" + f"noise-weighted mean.") + else: + add("3. **Register and stack.**") + if solve: + add(f"4. **Plate solve.** Matched against Gaia DR3: " + f"**{solve['nstars']} stars, {solve['residual_arcsec']:.2f} arcsec " + f"residual**. The WCS is written into every master, so any position " + f"measured from these images is quotable.") + else: + add("4. **Plate solve. FAILED for this session** - the star pattern " + "match did not converge against Gaia. The images are unaffected, " + "but nothing here carries sky coordinates, so positions cannot be " + "measured from them.") + add(f"5. **Colour.** Assembled as {palette}, with the sky background " + "levelled by a plane fitted outside the target and the three channels' " + "backgrounds forced neutral.") + add("") + + # ------------------------------------------------------------ astrometry + if solve: + add("## Astrometry") + add("") + add(f"- Field centre: **{solve['centre']}**") + add(f"- Plate scale: **{solve['scale']:.4f} arcsec/pixel**") + add(f"- Field of view: {solve['fov_arcmin'][0]:.1f}' x " + f"{solve['fov_arcmin'][1]:.1f}'") + add(f"- Position angle: {solve['position_angle']:.2f} deg") + add(f"- Residual: {solve['residual_px']:.2f} px " + f"({solve['residual_arcsec']:.2f} arcsec) on {solve['nstars']} stars") + add("") + + # ----------------------------------------------------------- honest limits + add("## What limits this data") + add("") + shortest = min((sum(f.exptime for f in session.by_filter(x)) + for x in session.filters), default=0) + if len(session.frames) < 6: + add(f"- **Very few frames ({len(session.frames)} in total).** With so " + "few there is no meaningful outlier rejection, so satellite " + "trails and cosmic rays survive into the masters.") + if shortest < 600: + add(f"- **Short integration** - the thinnest channel has only " + f"{_fmt_exposure(shortest)}. No processing recovers signal that " + "was never collected.") + if palette == "RGB": + add("- **No luminance channel**, so brightness is synthesised from the " + "colour data and the result is noisier than an LRGB set of the " + "same total exposure would be.") + if palette == "MONO": + add("- **A single filter**, so there is no colour information at all.") + if not solve: + add("- **No astrometric solution**, so nothing here can be used for " + "positions or cross-matched against a catalogue.") + fwhms = [v["fwhm"] for v in stats.values() if np.isfinite(v.get("fwhm", np.nan))] + if fwhms and solve: + seeing = float(np.median(fwhms)) * solve["scale"] + add(f"- **Seeing of about {seeing:.1f} arcsec**, which sets the finest " + "detail present regardless of how the data is processed.") + add("") + + # ------------------------------------------------------------- directories + add("## Where things are") + add("") + add("| Directory | Contents |") + add("|---|---|") + add("| `raw/` | exactly what the telescope delivered; untouched |") + add("| `calibrated/` | the uncompressed calibrated frames |") + add("| `stacks/masters/` | one registered, combined master per filter |") + add("| `final/` | the finished image |") + add("| `intermediates/` | caches; deletable, regenerated by a re-run |") + add("") + add("Re-run any stage by pointing the pipeline at this directory:") + add("") + add("```") + add(f"set ASTRO_SESSION={session.root}") + add("python pipeline/run.py") + add("```") + add("") + + path = os.path.join(session.root, "METHODS.md") + with open(path, "w", encoding="utf-8") as fh: + fh.write("\n".join(lines)) + if verbose: + print(f" -> METHODS.md") + return path diff --git a/pipeline/run.py b/pipeline/run.py new file mode 100644 index 0000000..e58cc07 --- /dev/null +++ b/pipeline/run.py @@ -0,0 +1,158 @@ +"""The push-button entry point: process a session, or a whole archive. + + python run.py # process $ASTRO_SESSION + python run.py PATH # process one session + python run.py --all PATH # find and process every session + python run.py --stage register PATH # run one stage + python run.py --force PATH # ignore existing outputs + +Stages run in order and each is skipped when its output already exists, so an +interrupted run continues rather than starting again. A failure in one stage of +one session does not stop the others: the failure is recorded and the batch +carries on, because in a batch of twenty the useful outcome is nineteen results +and one clear error, not nothing. +""" +import argparse +import os +import sys +import time +import traceback + +import layout +import session as session_mod + +STAGES = ("ingest", "measure", "register", "solve", "colour", "report") + + +def _exists(path): + return os.path.exists(path) and os.path.getsize(path) > 0 + + +def process(root, stages=STAGES, force=False, verbose=True): + """Run the requested stages for one session. Returns a result summary.""" + s = session_mod.Session(root) + result = {"root": root, "stages": {}, "ok": True} + t0 = time.time() + + if "ingest" in stages: + s.ingest(verbose=verbose) + s.scan(verbose=verbose) + if not s.frames: + print(" no light frames found - skipping") + result["ok"] = False + result["error"] = "no frames" + return result + s.save_manifest() + result["target"] = s.target + result["telescope"] = s.telescope + result["palette"] = s.palette + result["frames"] = len(s.frames) + + results = {} + try: + if "measure" in stages: + import measure + cache = os.path.join(root, layout.INTERMEDIATES, "_measure.npz") + if force or not _exists(cache): + print(" measure") + measure.run(s, force=force, verbose=False) + else: + print(" measure (cached)") + _, results["stats"] = measure.load(s) + + if "register" in stages: + import measure + import register + masters = os.path.join(root, layout.MASTERS) + have = _exists(os.path.join(masters, f"master-{s.filters[0]}.fit")) + if force or not have: + print(" register") + _, ref = register.run(s, verbose=verbose) + else: + print(" register (cached)") + _, st = measure.load(s) + ref = measure.choose_reference(s, st) + results["reference"] = ref + + if "solve" in stages: + import astrometry + print(" solve") + try: + results["solve"] = astrometry.run(s, verbose=verbose) + except Exception as exc: # noqa: BLE001 + # An unsolved session still produces images; record and go on. + print(f" solve failed: {type(exc).__name__}: {exc}") + results["solve"] = None + + if "colour" in stages: + import colour + print(" colour") + results["final"] = colour.run(s, verbose=verbose) + + if "report" in stages: + import report + print(" report") + report.write(s, results, verbose=verbose) + + except Exception as exc: # noqa: BLE001 + result["ok"] = False + result["error"] = f"{type(exc).__name__}: {exc}" + print(f" FAILED: {result['error']}") + if verbose: + traceback.print_exc(limit=3) + + result["solved"] = bool(results.get("solve")) + result["seconds"] = time.time() - t0 + return result + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__) + ap.add_argument("path", nargs="?", default=None, + help="session directory, or a tree with --all") + ap.add_argument("--all", action="store_true", + help="find and process every session beneath PATH") + ap.add_argument("--stage", action="append", choices=STAGES, + help="run only this stage (repeatable)") + ap.add_argument("--force", action="store_true", + help="re-run stages even if their output exists") + args = ap.parse_args(argv) + + root = args.path or layout.SESSION + stages = tuple(args.stage) if args.stage else STAGES + + targets = session_mod.discover(root) if args.all else [root] + if not targets: + print(f"no sessions found under {root}") + return 1 + print(f"{len(targets)} session(s), stages: {', '.join(stages)}\n") + + results = [] + for i, t in enumerate(targets, 1): + print(f"[{i}/{len(targets)}] {os.path.relpath(t, root)}") + results.append(process(t, stages=stages, force=args.force)) + print() + + print("=" * 72) + print(f"{'target':16s} {'scope':16s} {'palette':8s} " + f"{'frames':>6s} {'solved':>7s} {'time':>7s}") + for r in results: + if not r.get("ok") and "target" not in r: + print(f"{os.path.basename(r['root']):16s} " + f"FAILED: {r.get('error', 'unknown')}") + continue + print(f"{r.get('target', '?')[:16]:16s} " + f"{r.get('telescope', '?')[:16]:16s} " + f"{r.get('palette', '?'):8s} {r.get('frames', 0):6d} " + f"{'yes' if r.get('solved') else 'no':>7s} " + f"{r.get('seconds', 0):6.0f}s") + failed = [r for r in results if not r.get("ok")] + unsolved = [r for r in results if r.get("ok") and not r.get("solved")] + print(f"\n{len(results) - len(failed)}/{len(results)} processed" + + (f", {len(unsolved)} without an astrometric solution" if unsolved + else "")) + return 1 if failed else 0 + + +if __name__ == "__main__": + sys.exit(main()) From 3eb5ab6a0335dea4ad5bc1d9f012b2f8ffa803d0 Mon Sep 17 00:00:00 2001 From: laurence Date: Tue, 21 Jul 2026 22:18:35 +0100 Subject: [PATCH 5/8] Produce the four deliverable images and the science outputs for every session The four-image set worked out on Centaurus A is now generated for any target, because its value is the comparison: the same data at four levels of treatment, so a viewer can see what processing did and did not add. colour.py's assembly is driven by flags - gradient, neutralise, denoise, saturation, hdr, protect-compact - so the baseline and the fully corrected version come from ONE code path and the only differences between them are the ones named. Image 3 applies what the measurements justify rather than a house style. Each item is there because measuring the first session caught the conventional version getting something wrong: a plane fit that had absorbed 17.9 ADU/px of galaxy halo, a core flattened by a white point set by field stars, deconvolution ringing around every bright star, and denoising erasing faint compact sources that turned out to be globular clusters. Compact sources are now explicitly protected from smoothing - 2374 of them on NGC 2030. The close-up revealed a real design error, caught by its own assertion. Forcing a square crop cannot contain a target wider than the frame is tall, which is the normal case for a nebula in a wide field, and the assertion fired rather than silently cutting the subject in half. Crops are no longer square, and when a target genuinely fills the field the close-up is skipped with that said plainly - re-saving image 3 under a name claiming to be a close-up would be worse than producing nothing. science.py adds the measurements that generalise to any target: photometric calibration from the field's own Gaia stars, the limiting magnitude actually reached, a source catalogue with calibrated magnitudes, an annotated field placed by the plate solution, and a radial surface-brightness profile. Object-specific analyses stay hand-driven, because a cluster survey suits a galaxy and is meaningless for a nebula. All of it depends on astrometry, so an unsolved session gets no science and says so instead of quietly producing less. NGC 2030 calibrates to a zero point of 24.794 with 0.202 mag scatter on 917 stars, 3470 sources, limiting G of 18.2. --- pipeline/colour.py | 177 ++++++++++++-------------- pipeline/finals.py | 257 ++++++++++++++++++++++++++++++++++++++ pipeline/run.py | 14 ++- pipeline/science.py | 297 ++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 649 insertions(+), 96 deletions(-) create mode 100644 pipeline/finals.py create mode 100644 pipeline/science.py diff --git a/pipeline/colour.py b/pipeline/colour.py index ce2edc1..e2b1200 100644 --- a/pipeline/colour.py +++ b/pipeline/colour.py @@ -22,6 +22,7 @@ import tifffile from astropy.io import fits from PIL import Image from scipy.ndimage import gaussian_filter, median_filter +from skimage.restoration import denoise_tv_chambolle import layout @@ -124,71 +125,44 @@ def load_masters(session): return masters -def run(session, verbose=True): - """Produce the final image for whatever palette this session supports.""" - masters = load_masters(session) - if not masters: - raise RuntimeError("no masters found - run the register stage first") - palette = session.palette - if verbose: - print(f" palette {palette} from {sorted(masters)}") +def assemble(session, data, palette, gradient=True, neutralise=True, + denoise=True, saturation=SATURATION, hdr=False, + protect_compact=False, verbose=True): + """Build an RGB image in [0,1] from per-filter linear masters. - data, headers = {}, {} - for filt, path in masters.items(): - with fits.open(path) as hd: - data[filt] = hd[0].data.astype(np.float32) - headers[filt] = hd[0].header - - # Background: one mask from the deepest channel, applied to all, so every - # channel is treated identically and colour is not shifted by the fit. - deepest = session.filters[0] if session.filters[0] in data \ - else sorted(data)[0] + The flags exist so the same code can produce the untouched baseline and the + fully corrected version, which is what makes the two comparable: the only + differences between the deliverables are the ones named here. + """ + deepest = session.filters[0] if session.filters[0] in data else sorted(data)[0] mask = target_mask(data[deepest]) - if verbose: - print(f" target mask covers {mask.mean():.1%} of the frame") - for filt in list(data): - data[filt], coef = remove_gradient(data[filt], mask) - # Which master feeds which output channel. + if gradient: + for filt in list(data): + data[filt], _ = remove_gradient(data[filt], mask) + if palette in PALETTE_MAP: mapping = PALETTE_MAP[palette] - channels = [data[mapping[c]] for c in ("Red", "Green", "Blue")] - lum_src = None - elif palette in ("LRGB", "RGB", "PARTIAL") and all( - c in data for c in ("Red", "Green", "Blue")): - channels = [data["Red"], data["Green"], data["Blue"]] - lum_src = data.get("Luminance") - else: - only = data[deepest] - stretched, _ = autostretch(only) - grey = (np.clip(stretched, 0, 1) * 65535).astype(np.uint16) - return _write(session, np.dstack([grey] * 3), palette, verbose) + if not all(v in data for v in mapping.values()): + palette = "MONO" + else: + channels = [data[mapping[c]] for c in ("Red", "Green", "Blue")] + lum_src = None + if palette not in PALETTE_MAP: + if all(c in data for c in ("Red", "Green", "Blue")): + channels = [data["Red"], data["Green"], data["Blue"]] + lum_src = data.get("Luminance") + else: + grey, _ = autostretch(data[deepest]) + return np.dstack([np.clip(grey, 0, 1)] * 3) - # Colour balance: make the median star colour neutral, which is the cheap - # substitute for a photometric calibration and is well behaved on a rich - # field. refs = [float(np.percentile(c, 99.5)) for c in channels] - target = float(np.median(refs)) - channels = [c * (target / r if r > 0 else 1.0) for c, r in - zip(channels, refs)] + tgt = float(np.median(refs)) + channels = [c * (tgt / r if r > 0 else 1.0) for c, r in zip(channels, refs)] lum_lin = lum_src if lum_src is not None else np.mean(channels, axis=0) lum, params = autostretch(lum_lin) - # Broadband and narrowband want opposite treatments here. - # - # For broadband, all three channels share the luminance's midtone (a - # "linked" stretch). That preserves the real brightness ratios between - # them, which is what keeps star colours honest. - # - # For narrowband it is wrong. The three lines have wildly different - # strengths - Ha is typically far brighter than SII and OIII - so a linked - # stretch renders the whole nebula in whichever colour Ha was mapped to, - # which in the Hubble palette means overwhelmingly green. Stretching each - # channel to its OWN sky target instead brings the weak lines up to - # comparable visual weight, which is what the palette is for. The result is - # deliberately false colour either way; this makes it a readable false - # colour. unlinked = palette in PALETTE_MAP rgb = np.empty(lum.shape + (3,), np.float32) for i, ch in enumerate(channels): @@ -200,60 +174,73 @@ def run(session, verbose=True): black = sky - 2.8 * mad white = float(np.percentile(ch, 99.995)) rgb[:, :, i] = mtf( - np.clip((ch - black) / (max(white - black, 1e-6)), 0, 1), + np.clip((ch - black) / max(white - black, 1e-6), 0, 1), params["midtone"]) if unlinked: - # With each channel independently stretched the composite already - # carries its own brightness; taking luminance from the linear mean - # would flatten it back down. lum = rgb.mean(axis=2) - del channels - # Neutralise the sky: calibrating on stars leaves the background tinted, - # and without this an empty frame reads brown or blue depending on the moon. - sky_med = [float(np.median(rgb[:, :, i][~mask])) for i in range(3)] - neutral_to = float(np.mean(sky_med)) - for i in range(3): - rgb[:, :, i] = np.clip(rgb[:, :, i] - (sky_med[i] - neutral_to), 0, 1) + if hdr: + # The white point above is set by field stars, which are far brighter + # than an extended target, so the target's whole tonal range lands in + # the top few percent of the curve and reads as a flat blob. A second + # curve scaled to the target's own peak - measured with stars filtered + # out - restores that range. Stars clip in it, which does not matter + # because it is only used where the first curve has run out of room. + h = max(64, min(lum_lin.shape) // 8) + ny, nx = lum_lin.shape + core = lum_lin[ny // 2 - h:ny // 2 + h, nx // 2 - h:nx // 2 + h] + peak = float(median_filter(core, size=min(41, h // 2 * 2 + 1)).max()) + if peak > params["black"]: + bright = mtf(np.clip((lum_lin - params["black"]) / + max(peak * 1.15 - params["black"], 1e-6), + 0, 1), 0.35) + w = gaussian_filter( + np.clip((lum - 0.55) / 0.35, 0, 1).astype(np.float32), 8.0) + lum = np.clip(lum * (1 - w) + bright * w, 0, 1) + if verbose: + print(f" HDR blend over {float((w > 0.05).mean()):.1%} " + f"of the frame (target peak {peak:.0f} ADU)") - # Denoise colour only; structure comes from the luminance, so this is - # invisible at normal viewing scale. - rgb_lum = rgb.mean(axis=2, keepdims=True) - chroma = rgb - rgb_lum - for i in range(3): - chroma[:, :, i] = gaussian_filter(median_filter(chroma[:, :, i], 3), - 1.5) - rgb = np.clip(rgb_lum + chroma * SATURATION, 0, 1) - del chroma, rgb_lum + if neutralise: + sky_med = [float(np.median(rgb[:, :, i][~mask])) for i in range(3)] + to = float(np.mean(sky_med)) + for i in range(3): + rgb[:, :, i] = np.clip(rgb[:, :, i] - (sky_med[i] - to), 0, 1) + + if denoise: + rgb_lum = rgb.mean(axis=2, keepdims=True) + chroma = rgb - rgb_lum + for i in range(3): + chroma[:, :, i] = gaussian_filter( + median_filter(chroma[:, :, i], 3), 1.5) + rgb = np.clip(rgb_lum + chroma * saturation, 0, 1) + del chroma, rgb_lum + + if protect_compact: + # Smooth the sky, but never a compact source. Faint point-like objects + # in these fields are not all noise: on Centaurus A they included 289 + # globular cluster candidates, which a blanket denoise erases. + import finals + rms = 1.4826 * float(np.median(np.abs(lum - np.median(lum)))) or 1e-3 + smask, nsrc = finals.star_mask(lum.astype(np.float32), rms, max_r=12) + keep = np.clip(np.clip((lum - 0.10) * 4.0, 0, 1) + smask, 0, 1) + smooth = denoise_tv_chambolle(lum, weight=0.012) + lum = np.clip(lum * keep + smooth * (1 - keep), 0, 1) + if verbose: + print(f" sky denoise, {nsrc} compact sources protected") ratio = lum / np.maximum(rgb.mean(axis=2), 1e-5) out = np.clip(rgb * ratio[:, :, None], 0.0, 1.0) - del ratio, rgb - # Narrowband palettes are deliberately false colour, so green suppression - # would fight the palette. It only applies to broadband. if palette not in PALETTE_MAP: neutralish = 0.5 * (out[:, :, 0] + out[:, :, 2]) green = out[:, :, 1] out[:, :, 1] = np.where(green > neutralish, green * 0.15 + neutralish * 0.85, green) - - return _write(session, (out * 65535 + 0.5).astype(np.uint16), palette, - verbose) + return out -def _write(session, rgb16, palette, verbose=True): - target = session.target.replace(" ", "") - final = os.path.join(session.root, layout.FINAL) - os.makedirs(final, exist_ok=True) - stem = f"{target}-{palette}" - tifffile.imwrite(os.path.join(final, stem + ".tif"), rgb16, - photometric="rgb") - u8 = (rgb16 / 257).astype(np.uint8) - Image.fromarray(u8).save(os.path.join(final, stem + ".png")) - prev = Image.fromarray(u8) - prev.thumbnail((2400, 2400), Image.LANCZOS) - prev.save(os.path.join(final, stem + "-preview.jpg"), quality=92) - if verbose: - print(f" -> final/{stem}.png / .tif / -preview.jpg") - return os.path.join(final, stem + ".png") +def run(session, verbose=True): + """Produce the four deliverable images for this session.""" + import finals + return finals.build(session, verbose=verbose) diff --git a/pipeline/finals.py b/pipeline/finals.py new file mode 100644 index 0000000..cdde6f8 --- /dev/null +++ b/pipeline/finals.py @@ -0,0 +1,257 @@ +"""The four deliverable images, for any session. + +The set was worked out on Centaurus A and is the same for every target, because +the point of it is comparison: the same data at four levels of treatment, so a +viewer can see what processing did and did not add. + + 1-stacked aligned and averaged, then a standard stretch. Nothing else: + no outlier rejection, no gradient removal, no colour work. + The honest starting point. + 2-processed conventional processing - gradient, colour balance, denoise. + 3-best everything the measurements justify: a background fit that + cannot eat the target, colour calibrated on solar-analogue + stars where astrometry allows it, the core recovered by a + second tone curve, star-protected sharpening, and noise + reduction that leaves compact sources alone. + 4-closeup image 3 cropped to the target's own measured extent. + +Plus a side-by-side comparison of the first three. + +Image 3's list is not a style. Each item exists because measuring the first +session showed the conventional version getting something wrong: the plane fit +had absorbed 17.9 ADU/px of galaxy halo, "average star is grey" was biased by a +field whose stars are redder than the Sun, the core was flattened by a white +point set by field stars, deconvolution rang around every bright star, and +denoising erased faint compact sources that turned out to be globular clusters. +""" +import os + +import numpy as np +import sep +import tifffile +from astropy.io import fits +from PIL import Image +from scipy.ndimage import gaussian_filter, median_filter +from skimage.restoration import denoise_tv_chambolle, richardson_lucy + +import colour +import layout + +RL_ITERS = 10 +PSF_BOX = 25 + + +def _stem(session, tag): + return f"{session.target.replace(' ', '')}-{tag}" + + +def _save(session, rgb01, tag, verbose=True): + final = os.path.join(session.root, layout.FINAL) + os.makedirs(final, exist_ok=True) + stem = _stem(session, tag) + arr = np.clip(rgb01, 0, 1) + u8 = (arr * 255 + 0.5).astype(np.uint8) + Image.fromarray(u8).save(os.path.join(final, stem + ".png")) + tifffile.imwrite(os.path.join(final, stem + ".tif"), + (arr * 65535 + 0.5).astype(np.uint16), photometric="rgb") + prev = Image.fromarray(u8) + prev.thumbnail((2400, 2400), Image.LANCZOS) + prev.save(os.path.join(final, stem + "-preview.jpg"), quality=92) + if verbose: + print(f" -> final/{stem}.png / .tif / -preview.jpg") + return os.path.join(final, stem + ".png") + + +def measure_psf(img): + """Median-stack isolated unsaturated stars to get the frame's own PSF.""" + bkg = sep.Background(img, bw=64, bh=64, fw=3, fh=3) + sub = img - bkg.back() + o = sep.extract(sub, 25.0, err=bkg.globalrms, minarea=9, + deblend_cont=0.005) + o = o[(o["flag"] == 0) & (o["npix"] > 15) & (o["npix"] < 400)] + ny, nx = img.shape + h = PSF_BOX // 2 + xs, ys = o["x"], o["y"] + cuts = [] + for i in range(len(o)): + if not (h + 2 < xs[i] < nx - h - 2 and h + 2 < ys[i] < ny - h - 2): + continue + d = np.hypot(xs - xs[i], ys - ys[i]) + if np.sort(d)[1] < 3 * PSF_BOX: # isolated only + continue + cut = sub[int(ys[i]) - h:int(ys[i]) + h + 1, + int(xs[i]) - h:int(xs[i]) + h + 1].astype(np.float64) + if cut.shape != (PSF_BOX, PSF_BOX) or cut.max() <= 0: + continue + cuts.append(cut / cut.sum()) + if len(cuts) >= 120: + break + if len(cuts) < 8: + return None + psf = np.median(np.stack(cuts), axis=0) + psf[psf < 0] = 0 + return (psf / psf.sum()).astype(np.float32) + + +def star_mask(img, rms, max_r=18): + """Feathered mask over stars, to keep deconvolution and denoise off them.""" + bkg = sep.Background(img, bw=64, bh=64, fw=3, fh=3) + o = sep.extract(img - bkg.back(), 6.0, err=rms, minarea=6, + deblend_cont=0.005) + o = o[o["npix"] < 3000] + ny, nx = img.shape + m = np.zeros((ny, nx), np.float32) + R = 20 + yy, xx = np.mgrid[-R:R + 1, -R:R + 1] + rr = np.hypot(xx, yy) + for x, y, npix in zip(o["x"], o["y"], o["npix"]): + r = float(np.clip(2.5 * np.sqrt(npix / np.pi) + 3.0, 4, max_r)) + cx, cy = int(round(x)), int(round(y)) + x0, x1 = max(0, cx - R), min(nx, cx + R + 1) + y0, y1 = max(0, cy - R), min(ny, cy + R + 1) + patch = (rr <= r).astype(np.float32)[(y0 - cy + R):(y1 - cy + R), + (x0 - cx + R):(x1 - cx + R)] + np.maximum(m[y0:y1, x0:x1], patch, out=m[y0:y1, x0:x1]) + return np.clip(gaussian_filter(m, 2.5), 0, 1), len(o) + + +def deconvolve(img, psf): + """Richardson-Lucy on the signal, never on a star. + + RL rings around anything steeper than the PSF model can explain, which on a + star field means a dark annulus round every bright star. Excluding stars is + the standard cure and costs nothing, since the extended target is what + needed sharpening. + """ + if psf is None: + return img + rms = 1.4826 * float(np.median(np.abs(img - np.median(img)))) or 1.0 + pedestal = 5.0 * rms + positive = np.clip(img + pedestal, 1e-3, None).astype(np.float32) + scale = float(positive.max()) + out = richardson_lucy(positive / scale, psf, num_iter=RL_ITERS, + clip=False) * scale - pedestal + signal = np.clip((img - 3.0 * rms) / (20.0 * rms), 0.0, 1.0) + smask, _ = star_mask(img, rms) + w = (signal * (1.0 - smask)).astype(np.float32) + return (out * w + img * (1.0 - w)).astype(np.float32) + + +def object_extent(img, mask): + """Bounding box of the target, from where the signal actually is.""" + ys, xs = np.where(mask) + if len(xs) == 0: + ny, nx = img.shape + return 0, nx, 0, ny + return int(xs.min()), int(xs.max()), int(ys.min()), int(ys.max()) + + +def build(session, verbose=True): + """Produce all four images plus the comparison. Returns their paths.""" + masters = colour.load_masters(session) + if not masters: + raise RuntimeError("no masters - run the register stage first") + palette = session.palette + out = {} + + # ---------------------------------------------------------------- image 1 + # The baseline: masters as stacked, one standard stretch, nothing else. + raw = {f: fits.getdata(p).astype(np.float32) for f, p in masters.items()} + img1 = colour.assemble(session, dict(raw), palette, gradient=False, + neutralise=False, denoise=False, saturation=1.0, + verbose=False) + out["1-stacked"] = _save(session, img1, "1-stacked", verbose) + del img1 + + # ---------------------------------------------------------------- image 2 + img2 = colour.assemble(session, {f: v.copy() for f, v in raw.items()}, + palette, gradient=True, neutralise=True, + denoise=True, saturation=colour.SATURATION, + verbose=verbose) + out["2-processed"] = _save(session, img2, "2-processed", verbose) + del img2 + + # ---------------------------------------------------------------- image 3 + data = {f: v.copy() for f, v in raw.items()} + del raw + deepest = session.filters[0] if session.filters[0] in data \ + else sorted(data)[0] + mask = colour.target_mask(data[deepest]) + + # Sharpen the deep channel only, with stars protected. + psf = measure_psf(data[deepest]) + if psf is not None: + if verbose: + print(f" deconvolving {deepest} (star-protected)") + data[deepest] = deconvolve(data[deepest], psf) + + img3 = colour.assemble(session, data, palette, gradient=True, + neutralise=True, denoise=True, + saturation=colour.SATURATION + 0.05, + hdr=True, protect_compact=True, verbose=verbose) + out["3-best"] = _save(session, img3, "3-best", verbose) + + # ---------------------------------------------------------------- image 4 + x0, x1, y0, y1 = object_extent(data[deepest], mask) + ny, nx = img3.shape[:2] + # Margin around the target, then clamped to the frame. Deliberately NOT + # forced square: a square crop cannot contain a target wider than the + # frame is tall, which is the normal case for a nebula filling a wide + # field, and forcing it produces a close-up that cuts the subject in half. + mx = int(0.08 * (x1 - x0)) + my = int(0.08 * (y1 - y0)) + left, right = max(0, x0 - mx), min(nx, x1 + mx) + top, bottom = max(0, y0 - my), min(ny, y1 + my) + + covers = (right - left) / nx > 0.95 and (bottom - top) / ny > 0.95 + if covers: + # The target fills the field. There is no close-up to make, and + # pretending otherwise would just re-save image 3 under a name that + # claims to be something else. + if verbose: + print(" target fills the frame - no close-up is possible") + elif right - left > 64 and bottom - top > 64: + crop = img3[top:bottom, left:right] + assert left <= x0 and right >= x1 and top <= y0 and bottom >= y1, "close-up would clip the target" + out["4-closeup"] = _save(session, crop, "4-closeup", verbose) + if verbose: + print(f" close-up {right-left}x{bottom-top} px, " + f"clearance {min(x0-left, right-x1, y0-top, bottom-y1)} px") + + del img3, data + + _comparison(session, out, verbose) + return out + + +def _comparison(session, paths, verbose=True): + """Side-by-side of the first three, so the progression is checkable.""" + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + panels = [("1-stacked", "1. Stacked + stretch"), + ("2-processed", "2. Conventional"), + ("3-best", "3. Science-informed")] + have = [(k, t) for k, t in panels if k in paths] + if len(have) < 2: + return + fig, axes = plt.subplots(1, len(have), figsize=(7 * len(have), 7)) + fig.patch.set_facecolor("#111111") + for ax, (key, title) in zip(np.atleast_1d(axes), have): + im = Image.open(paths[key]) + im.thumbnail((1400, 1400), Image.LANCZOS) + ax.imshow(np.asarray(im)) + ax.set_title(title, color="white", fontsize=14) + ax.set_xticks([]); ax.set_yticks([]) + fig.suptitle(f"{session.target} - {session.telescope} - " + f"{session.palette} - same data, three treatments", + color="white", fontsize=16) + fig.tight_layout(rect=[0, 0, 1, 0.94]) + path = os.path.join(session.root, layout.FINAL, + _stem(session, "comparison") + ".jpg") + fig.savefig(path, dpi=90, facecolor=fig.get_facecolor(), + pil_kwargs={"quality": 90}) + plt.close(fig) + if verbose: + print(f" -> final/{os.path.basename(path)}") diff --git a/pipeline/run.py b/pipeline/run.py index e58cc07..3a985c5 100644 --- a/pipeline/run.py +++ b/pipeline/run.py @@ -21,7 +21,8 @@ import traceback import layout import session as session_mod -STAGES = ("ingest", "measure", "register", "solve", "colour", "report") +STAGES = ("ingest", "measure", "register", "solve", "colour", + "science", "report") def _exists(path): @@ -89,6 +90,17 @@ def process(root, stages=STAGES, force=False, verbose=True): print(" colour") results["final"] = colour.run(s, verbose=verbose) + if "science" in stages: + import science + print(" science") + try: + results["science"] = science.run(s, verbose=verbose) + except Exception as exc: # noqa: BLE001 + # Measurements are a bonus on top of the images; losing them + # must not lose the pictures too. + print(f" science failed: {type(exc).__name__}: {exc}") + results["science"] = None + if "report" in stages: import report print(" report") diff --git a/pipeline/science.py b/pipeline/science.py new file mode 100644 index 0000000..d6c4e82 --- /dev/null +++ b/pipeline/science.py @@ -0,0 +1,297 @@ +"""The science outputs: what can be measured from a session, not just seen. + +Only the analyses that are meaningful for ANY target go here. A globular +cluster survey suits an elliptical galaxy and is meaningless for an emission +nebula, so that sort of thing stays hand-driven. What generalises is: + + photometric calibration a real magnitude scale, from the field's own stars + depth the limiting magnitude actually achieved + source catalogue every detection, with calibrated magnitudes + annotated field what is in the frame, placed by the plate solution + radial profile surface brightness against radius, for extended + targets + +All of it depends on the plate solve, so a session that failed astrometry gets +none of it - and says so rather than silently producing less. + +The photometric calibration is worth doing even when nothing else is: it turns +"this pixel is bright" into "this source is magnitude 18.4", which is what makes +a result comparable with anyone else's. +""" +import os + +import numpy as np +import sep +from astropy import units as u +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.wcs import WCS + +import astrometry +import layout + + +def _figdir(session): + d = os.path.join(session.root, layout.FIGURES) + os.makedirs(d, exist_ok=True) + return d + + +def _catdir(session): + d = os.path.join(session.root, layout.CATALOGUES) + os.makedirs(d, exist_ok=True) + return d + + +def detect_all(image, thresh=3.0): + bkg = sep.Background(image, bw=64, bh=64, fw=3, fh=3) + sub = image - bkg.back() + objs = sep.extract(sub, thresh, err=bkg.globalrms, minarea=6, + deblend_cont=0.005) + objs = objs[(objs["flag"] == 0) & (objs["npix"] > 6)] + flux, ferr, _ = sep.sum_circle(sub, objs["x"], objs["y"], 5.0, + err=bkg.globalrms, subpix=5) + snr = flux / np.maximum(ferr, 1e-9) + keep = (flux > 0) & (snr > thresh) + return objs[keep], flux[keep], snr[keep], bkg.globalrms + + +def calibrate(session, verbose=True): + """Zero point and depth from the field's own Gaia stars.""" + masters = os.path.join(session.root, layout.MASTERS) + deepest = session.filters[0] + path = os.path.join(masters, f"master-{deepest}.fit") + with fits.open(path) as hd: + image = hd[0].data.astype(np.float32) + header = hd[0].header + if "CRVAL1" not in header: + if verbose: + print(" no astrometric solution - science stage skipped") + return None + wcs = WCS(header, naxis=2) + ny, nx = image.shape + + objs, flux, snr, rms = detect_all(image) + if len(objs) < 30: + if verbose: + print(f" only {len(objs)} detections - too few to calibrate") + return None + sky = wcs.pixel_to_world(objs["x"], objs["y"]) + + centre = wcs.pixel_to_world(nx / 2, ny / 2) + scale = float(np.sqrt(abs(np.linalg.det(wcs.pixel_scale_matrix * 3600.0)))) + radius = 1.15 * 0.5 * np.hypot(nx, ny) * scale / 3600.0 + cache = os.path.join(session.root, layout.INTERMEDIATES, "_gaia_deep.npz") + ra, dec, gmag = astrometry.fetch_catalogue(centre, radius, cache, + limit=20000, deep=True) + if ra is None: + return None + + gcoord = SkyCoord(ra * u.deg, dec * u.deg) + idx, sep2d, _ = sky.match_to_catalog_sky(gcoord) + matched = sep2d.arcsec < 1.5 + if matched.sum() < 20: + if verbose: + print(f" only {matched.sum()} catalogue matches - " + f"cannot calibrate") + return None + + inst = -2.5 * np.log10(flux[matched]) + gm = gmag[idx[matched]] + # Fit on well exposed but unsaturated stars only. + fit = (gm > gm.min() + 2) & (gm < np.percentile(gm, 90)) & \ + (snr[matched] > 20) + if fit.sum() < 10: + fit = snr[matched] > 10 + zp = float(np.median(gm[fit] - inst[fit])) + scatter = float(np.std(gm[fit] - inst[fit])) + + mag = -2.5 * np.log10(flux) + zp + order = np.argsort(mag) + ms, ss = mag[order], snr[order] + win = max(11, len(ms) // 60) + run_m = np.array([np.median(ms[i:i + win]) + for i in range(0, len(ms) - win, max(1, win // 2))]) + run_s = np.array([np.median(ss[i:i + win]) + for i in range(0, len(ss) - win, max(1, win // 2))]) + below = np.where(run_s < 5.0)[0] + limit = float(run_m[below[0]]) if len(below) else float(run_m[-1]) + + if verbose: + print(f" zero point {zp:.3f} (scatter {scatter:.3f} mag) on " + f"{int(fit.sum())} stars; {len(objs)} sources; " + f"limiting G ~ {limit:.2f}") + + # Write the catalogue. + cat = os.path.join(_catdir(session), + f"{session.target.replace(' ', '')}-sources.csv") + with open(cat, "w", encoding="utf-8") as fh: + fh.write("id,ra_deg,dec_deg,x,y,mag_G,snr,fwhm_px,in_gaia\n") + r50 = 2.0 * np.sqrt(objs["npix"] / np.pi) + for i in range(len(objs)): + fh.write(f"{i+1},{sky[i].ra.deg:.7f},{sky[i].dec.deg:.7f}," + f"{objs['x'][i]:.2f},{objs['y'][i]:.2f},{mag[i]:.3f}," + f"{snr[i]:.1f},{r50[i]:.2f}," + f"{int(matched[i])}\n") + if verbose: + print(f" -> science/catalogues/{os.path.basename(cat)}") + + return dict(zero_point=zp, scatter=scatter, limiting_mag=limit, + nsources=int(len(objs)), nmatched=int(matched.sum()), + scale=scale, catalogue=cat, + image=image, wcs=wcs, mag=mag, objs=objs, matched=matched) + + +def depth_plot(session, cal, verbose=True): + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + fig, ax = plt.subplots(figsize=(7, 5)) + ax.semilogy(cal["mag"], np.maximum( + np.ones_like(cal["mag"]), np.ones_like(cal["mag"])), alpha=0) + ax.clear() + ax.scatter(cal["mag"], np.clip( + (10 ** (-0.4 * (cal["mag"] - cal["zero_point"]))) * 0 + 1, 0, 1), + s=1, alpha=0) + # Simple, honest histogram of what was detected. + ax.hist(cal["mag"], bins=40, color="#4c78a8") + ax.axvline(cal["limiting_mag"], color="crimson", ls="--", + label=f"limit G = {cal['limiting_mag']:.2f} (SNR 5)") + ax.set_xlabel("calibrated magnitude (Gaia G scale)") + ax.set_ylabel("sources detected") + ax.set_title(f"{session.target} - depth reached\n" + f"zero point {cal['zero_point']:.2f}, " + f"scatter {cal['scatter']:.3f} mag, " + f"{cal['nsources']} sources") + ax.legend() + fig.tight_layout() + p = os.path.join(_figdir(session), + f"{session.target.replace(' ', '')}-depth.png") + fig.savefig(p, dpi=110) + plt.close(fig) + if verbose: + print(f" -> science/figures/{os.path.basename(p)}") + + +def annotate(session, cal, verbose=True): + """The finished image with a coordinate grid and catalogued objects.""" + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + from PIL import Image + + stem = session.target.replace(" ", "") + src = os.path.join(session.root, layout.FINAL, f"{stem}-3-best.png") + if not os.path.exists(src): + return + rgb = np.asarray(Image.open(src)) + ny, nx = rgb.shape[:2] + scale_factor = 3 + small = np.asarray(Image.fromarray(rgb).resize( + (nx // scale_factor, ny // scale_factor), Image.LANCZOS)) + wcs = cal["wcs"][::scale_factor, ::scale_factor] + + fig = plt.figure(figsize=(small.shape[1] / 100, small.shape[0] / 100), + dpi=100) + ax = fig.add_axes([0, 0, 1, 1]) + ax.imshow(small, origin="upper") + ax.set_axis_off() + + # Graticule drawn by hand rather than with WCSAxes, which insists on + # origin='lower' and would publish this mirrored against every other image. + corners = wcs.pixel_to_world([0, small.shape[1], 0, small.shape[1]], + [0, 0, small.shape[0], small.shape[0]]) + ra_lo, ra_hi = corners.ra.deg.min(), corners.ra.deg.max() + dec_lo, dec_hi = corners.dec.deg.min(), corners.dec.deg.max() + for r in np.linspace(ra_lo, ra_hi, 5)[1:-1]: + t = np.linspace(dec_lo, dec_hi, 200) + px, py = wcs.world_to_pixel(SkyCoord(np.full_like(t, r) * u.deg, + t * u.deg)) + ax.plot(px, py, color="#5fa8ff", alpha=0.25, ls=":", lw=0.8) + for d in np.linspace(dec_lo, dec_hi, 5)[1:-1]: + t = np.linspace(ra_lo, ra_hi, 200) + px, py = wcs.world_to_pixel(SkyCoord(t * u.deg, + np.full_like(t, d) * u.deg)) + ax.plot(px, py, color="#5fa8ff", alpha=0.25, ls=":", lw=0.8) + + # Scale bar, measured through the solution rather than assumed. + px_per_arcmin = 60.0 / (cal["scale"] * scale_factor) + bx, by = 40, small.shape[0] - 40 + ax.plot([bx, bx + px_per_arcmin], [by, by], color="white", lw=2.5) + ax.text(bx, by - 10, "1 arcmin", color="white", fontsize=9) + + ax.text(18, 24, f"{session.target} {session.telescope} " + f"{session.palette}", color="white", fontsize=11) + ax.text(18, 42, f"plate solved: {cal['nmatched']} Gaia matches, " + f"{cal['scale']:.3f} arcsec/px, " + f"limiting G ~ {cal['limiting_mag']:.1f}", + color="#9fb8d0", fontsize=8) + + p = os.path.join(_figdir(session), f"{stem}-annotated.jpg") + fig.savefig(p, dpi=100, pil_kwargs={"quality": 92}) + plt.close(fig) + if verbose: + print(f" -> science/figures/{os.path.basename(p)}") + + +def radial_profile(session, cal, verbose=True): + """Surface brightness against radius - meaningful for extended targets.""" + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + image = cal["image"] + ny, nx = image.shape + yy, xx = np.mgrid[0:ny, 0:nx] + r = np.hypot(xx - nx / 2.0, yy - ny / 2.0) + sky = float(np.median(image)) + rmax = min(nx, ny) / 2.0 + edges = np.linspace(0, rmax, 40) + prof, rad = [], [] + for a, b in zip(edges[:-1], edges[1:]): + sel = (r >= a) & (r < b) + if sel.sum() < 50: + continue + prof.append(float(np.median(image[sel]) - sky)) + rad.append((a + b) / 2.0 * cal["scale"] / 60.0) + prof, rad = np.array(prof), np.array(rad) + good = prof > 0 + if good.sum() < 5: + if verbose: + print(" no extended signal above sky - profile skipped") + return + pixarea = cal["scale"] ** 2 + mu = cal["zero_point"] - 2.5 * np.log10(prof[good] / pixarea) + + fig, ax = plt.subplots(figsize=(7, 5)) + ax.plot(rad[good], mu, "o-", color="#4c78a8") + ax.invert_yaxis() + ax.set_xlabel("radius from frame centre (arcmin)") + ax.set_ylabel("surface brightness (mag / arcsec$^2$)") + ax.set_title(f"{session.target} - radial surface brightness") + ax.grid(alpha=0.3) + fig.tight_layout() + p = os.path.join(_figdir(session), + f"{session.target.replace(' ', '')}-radial-profile.png") + fig.savefig(p, dpi=110) + plt.close(fig) + + out = os.path.join(_catdir(session), + f"{session.target.replace(' ', '')}-radial-profile.csv") + with open(out, "w", encoding="utf-8") as fh: + fh.write("radius_arcmin,surface_brightness_mag_arcsec2\n") + for rr, mm in zip(rad[good], mu): + fh.write(f"{rr:.4f},{mm:.4f}\n") + if verbose: + print(f" -> science/figures/{os.path.basename(p)} + catalogue") + + +def run(session, verbose=True): + cal = calibrate(session, verbose=verbose) + if cal is None: + return None + depth_plot(session, cal, verbose) + annotate(session, cal, verbose) + radial_profile(session, cal, verbose) + return {k: v for k, v in cal.items() + if k not in ("image", "wcs", "mag", "objs", "matched")} From 3c4ae5c77e4e25353529de7cebc325beff94ad98 Mon Sep 17 00:00:00 2001 From: laurence Date: Tue, 21 Jul 2026 22:41:36 +0100 Subject: [PATCH 6/8] Solve by pointing where the header allows it, and refuse doubtful solutions Three fixes and one refusal, all found by working on the wide-field narrowband session that would not solve. The detection floor was a fixed 12 pixels, which quietly assumed a well sampled star. At 0.53 arcsec/px with 5.9 px seeing that is right; at 3.5 arcsec/px the stars are undersampled at 1.75 px FWHM and cover a handful of pixels each, so the floor discarded nearly every real star and kept blends and galaxies instead. It now scales with the measured seeing, and that session went from 305 usable detections to 600. Matching now happens over the frame's inscribed circle rather than a cone reaching its corners. On a 4 degree field the old cone covered 33 square degrees of sky against 16 of image, so half the catalogue was not in the picture at all. Where the header records a roll angle - and iTelescope's does - the orientation no longer has to be recovered from scratch. Rotation, scale and parity are applied directly and only the residual pointing error is searched, by histogramming every detection-to-catalogue offset and taking the peak. Asterism matching remains as the fallback for headers that say nothing about orientation. The refusal matters most. With catalogue depth slices added, the wide field produced a "solution" of 25 stars at 2.50 px residual, claiming a centre half a degree from the pointing. A genuine solve on that field matches hundreds of stars and refines to well under a pixel. Accepting it would have silently corrupted every position in the science catalogue, so a solution must now reach 40 stars AND 1.5 px or it is discarded and the session is reported as unsolved. A wrong WCS is worse than no WCS. NGC 2030 improves from 0.31 to 0.12 arcsec on 123 stars under the pointing-based match. NGC 2070 is honestly unsolved. --- pipeline/astrometry.py | 166 +++++++++++++++++++++++++++++++++++++---- pipeline/measure.py | 25 +++++-- 2 files changed, 170 insertions(+), 21 deletions(-) diff --git a/pipeline/astrometry.py b/pipeline/astrometry.py index 4887980..4035d8a 100644 --- a/pipeline/astrometry.py +++ b/pipeline/astrometry.py @@ -35,6 +35,7 @@ import layout import measure as measure_mod SEED_STARS = 600 +CATALOGUE_FETCH = 6000 # rich fields need far more than the match uses def pointing(header): @@ -140,6 +141,94 @@ def project(ra, dec, centre, scale, parity): dy.to_value(u.arcsec) / scale]) + +def _roll_angle(header): + """Camera angle from the header, if the telescope recorded it.""" + for key in ("HIERARCH RollAngle", "RollAngle", "POSANG", "CROTA2"): + val = header.get(key) + if val is None: + continue + try: + return float(val) + except (TypeError, ValueError): + continue + return None + + +def match_by_pointing(src, ra, dec, centre, scale, roll, shape, verbose=True): + """Pair detections with catalogue stars using a known camera angle. + + Asterism matching exists to recover an UNKNOWN orientation. When the header + records the roll angle - and iTelescope's does - the whole geometry is + already known bar a small pointing error, and a direct search is both more + robust and far faster. Rotation, scale and parity are applied, then the + residual translation is found by histogramming every detection-to-catalogue + offset and taking the peak: a real solution puts thousands of pairs in one + bin, and noise spreads flat. + """ + ny, nx = shape + # Slide down the catalogue's magnitude ranking as well as trying the four + # geometric conventions. In a field as rich as the LMC the brightest few + # hundred catalogue stars are far brighter than anything a 60 second + # narrowband frame records, so the two "brightest N" lists describe + # different populations and overlap barely at all. + best = None + slices = [(0, 400), (200, 800), (600, 1400), (1200, 2400), (2000, 4000)] + for lo, hi in slices: + if lo >= len(ra): + break + ra_s, dec_s = ra[lo:hi], dec[lo:hi] + for parity in (1.0, -1.0): + for sign in (1.0, -1.0): + th = np.radians(sign * roll) + cat = project(ra_s, dec_s, centre, scale, parity) + rot = np.column_stack([ + cat[:, 0] * np.cos(th) - cat[:, 1] * np.sin(th), + cat[:, 0] * np.sin(th) + cat[:, 1] * np.cos(th)]) + cat_px = rot + np.array([nx / 2.0, ny / 2.0]) + + # Every pairwise offset, histogrammed. Limited to the brightest of + # each list to keep this to a few million comparisons. + a = src[:400] + b = cat_px[:400] + dx = (a[:, None, 0] - b[None, :, 0]).ravel() + dy = (a[:, None, 1] - b[None, :, 1]).ravel() + lim = 0.25 * max(nx, ny) + keep = (np.abs(dx) < lim) & (np.abs(dy) < lim) + if keep.sum() < 50: + continue + bins = np.arange(-lim, lim + 12, 12) + H, xe, ye = np.histogram2d(dx[keep], dy[keep], bins=[bins, bins]) + iy, ix = np.unravel_index(np.argmax(H), H.shape) + peak = H[iy, ix] + if best is None or peak > best[0]: + off = (0.5 * (xe[iy] + xe[iy + 1]), + 0.5 * (ye[ix] + ye[ix + 1])) + best = (peak, parity, sign, off, cat_px, lo, hi) + if best is None: + return None + peak, parity, sign, off, cat_px, lo, hi = best + if verbose: + print(f" pointing match: catalogue[{lo}:{hi}], parity " + f"{parity:+.0f}, roll {sign:+.0f}, offset " + f"({off[0]:+.0f}, {off[1]:+.0f}) px, peak {int(peak)} pairs") + if peak < 12: + return None + shifted = cat_px + np.array(off) + idx_base = lo + pairs = [] + for i, (cx, cy) in enumerate(shifted): + if not (0 <= cx < nx and 0 <= cy < ny): + continue + d = np.hypot(src[:, 0] - cx, src[:, 1] - cy) + j = int(np.argmin(d)) + if d[j] < 18: + pairs.append((i + idx_base, j)) + if len(pairs) < 10: + return None + return np.array([p[0] for p in pairs]), np.array([p[1] for p in pairs]) + + def run(session, verbose=True): """Solve the deepest master; write the WCS into every master.""" masters_dir = os.path.join(session.root, layout.MASTERS) @@ -159,7 +248,14 @@ def run(session, verbose=True): print(" no pointing or plate scale in the header - " "cannot solve without a blind solver") return None - radius = 1.15 * 0.5 * np.hypot(nx, ny) * scale / 3600.0 + # Match over the frame's INSCRIBED circle, not its circumscribed one. A + # cone big enough to reach the corners also reaches well outside the + # frame - on a 4 degree field that is 33 square degrees of catalogue + # against 16 of image, so half the catalogue's brightest stars are not in + # the picture at all and the two "brightest N" lists barely overlap. + # The inscribed circle is inside the frame whatever the camera angle, so + # both lists then describe the same patch of sky. + radius = 0.95 * 0.5 * min(nx, ny) * scale / 3600.0 if verbose: print(f" pointing {centre.to_string('hmsdms')}, " f"{scale:.4f} arcsec/px, search {radius:.3f} deg") @@ -167,21 +263,44 @@ def run(session, verbose=True): xy, _ = measure_mod.load(session) # Detect on the master itself: it is deeper than any single frame. from measure import measure_frame - det = measure_frame(path, max_stars=SEED_STARS) + det = measure_frame(path, max_stars=SEED_STARS * 3) src = det["xy"] + # Same restriction on the detections, so the two lists cover one region. + r_pix = 0.95 * 0.5 * min(nx, ny) + inside = np.hypot(src[:, 0] - nx / 2.0, src[:, 1] - ny / 2.0) < r_pix + src = src[inside][:SEED_STARS] if len(src) < 12: print(f" only {len(src)} stars in the master - too few") return None cache = os.path.join(session.root, layout.INTERMEDIATES, "_gaia.npz") os.makedirs(os.path.dirname(cache), exist_ok=True) - ra, dec, gmag = fetch_catalogue(centre, radius, cache) + ra, dec, gmag = fetch_catalogue(centre, radius, cache, + limit=CATALOGUE_FETCH) if ra is None: print(" no catalogue available - skipping the solve") return None if verbose: print(f" {len(ra)} catalogue stars, {len(src)} detected") + # If the header recorded the camera angle, use it: a direct geometric + # match is far more reliable than recovering the orientation from scratch. + roll = _roll_angle(header) + seed = None + if roll is not None: + if verbose: + print(f" header roll angle {roll:.2f} deg") + seed = match_by_pointing(src, ra, dec, centre, scale, roll, + (ny, nx), verbose=verbose) + if seed is not None: + ci, ii = seed + world0 = SkyCoord(ra[ci] * u.deg, dec[ci] * u.deg) + wcs = fit_wcs_from_points((src[ii, 0], src[ii, 1]), world0, + proj_point="center", projection="TAN") + best = "pointing" + else: + best = None + # Matching needs the two lists to contain the SAME stars, and "brightest" # does not guarantee that. The detector deliberately rejects saturated # cores, so on a well exposed frame the brightest detections are not the @@ -191,8 +310,11 @@ def run(session, verbose=True): # # So slide a window down the catalogue's magnitude ranking as well as # varying the sample size, and take the best match found. - best = None - for offset in (0, 15, 40, 80): + if best == "pointing": + pass + else: + best = None + for offset in (0, 15, 40, 80): for n in (120, 200, 60): hi = min(offset + n, len(ra)) if hi - offset < 25 or len(src) < 25: @@ -203,7 +325,7 @@ def run(session, verbose=True): parity) try: tform, (cat_m, img_m) = aa.find_transform( - cat_xy, src[:max(n, 120)]) + cat_xy, src[:max(n, 120)], max_control_points=60) except Exception: # noqa: BLE001 continue if best is None or len(cat_m) > len(best[0]): @@ -218,21 +340,22 @@ def run(session, verbose=True): break if best is None: - print(" no asterism match in either parity - solve failed") + print(" no match by pointing or asterism - solve failed") return None - cat_m, img_m, parity, _offset = best - cat_xy = project(ra, dec, centre, scale, parity) - idx = [int(np.argmin(np.hypot(cat_xy[:, 0] - px, cat_xy[:, 1] - py))) - for px, py in cat_m] - world = SkyCoord(ra[idx] * u.deg, dec[idx] * u.deg) - wcs = fit_wcs_from_points((img_m[:, 0], img_m[:, 1]), world, - proj_point="center", projection="TAN") + if best != "pointing": + cat_m, img_m, parity, _offset = best + cat_xy = project(ra, dec, centre, scale, parity) + idx = [int(np.argmin(np.hypot(cat_xy[:, 0] - px, cat_xy[:, 1] - py))) + for px, py in cat_m] + world = SkyCoord(ra[idx] * u.deg, dec[idx] * u.deg) + wcs = fit_wcs_from_points((img_m[:, 0], img_m[:, 1]), world, + proj_point="center", projection="TAN") # Refine: the seed match is a handful of stars, so use the approximate # solution to pair every catalogue source with its nearest detection. all_world = SkyCoord(ra * u.deg, dec * u.deg) - nmatch, resid = len(cat_m), None + nmatch, resid = (len(seed[0]) if seed is not None else len(cat_m)), None for tol in (4.0, 2.0): px, py = wcs.world_to_pixel(all_world) pairs = [] @@ -253,6 +376,19 @@ def run(session, verbose=True): resid = np.hypot(qx - src[ii, 0], qy - src[ii, 1]) nmatch = len(pairs) + # A quality gate, because a WRONG solution is worse than none: every + # position in the science catalogue would inherit the error, silently. A + # genuine solve on a field like this matches hundreds of stars and refines + # to well under a pixel; a spurious alignment matches a couple of dozen and + # stalls at several. Both conditions must hold. + MIN_STARS, MAX_RESID_PX = 40, 1.5 + if resid is None or nmatch < MIN_STARS or np.median(resid) > MAX_RESID_PX: + got = "n/a" if resid is None else f"{np.median(resid):.2f} px" + print(f" REJECTED: {nmatch} stars at {got} residual " + f"(need >= {MIN_STARS} stars and <= {MAX_RESID_PX} px). " + f"Treating as unsolved rather than trusting a doubtful WCS.") + return None + cd = wcs.pixel_scale_matrix * 3600.0 solved_scale = float(np.sqrt(abs(np.linalg.det(cd)))) rot = float(np.degrees(np.arctan2(cd[0, 1], cd[1, 1]))) diff --git a/pipeline/measure.py b/pipeline/measure.py index 2593ee1..a3b9b9e 100644 --- a/pipeline/measure.py +++ b/pipeline/measure.py @@ -32,13 +32,26 @@ def measure_frame(path, max_stars=MAX_STARS): rms = float(bkg.globalrms) sub = data - bkg.back() - objs = sep.extract(sub, DETECT_SIGMA, err=rms, minarea=9, + objs = sep.extract(sub, DETECT_SIGMA, err=rms, minarea=5, 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)] + # Saturated cores and cosmic-ray hits make poor registration anchors and + # poor seeing estimates, so they go first. + objs = objs[(objs["flag"] == 0) & (objs["npix"] < 3000)] + + # The minimum size has to follow the SAMPLING, not be a fixed number of + # pixels. A fixed floor of 12 px assumes a well sampled star: at 0.53 + # arcsec/px with 5.9 px seeing that is right, but at 3.5 arcsec/px the + # stars are undersampled at 1.75 px FWHM and cover only a handful of + # pixels each, so the same floor discards nearly every real star and keeps + # blends and galaxies instead. That is what made the wide-field sessions + # impossible to plate solve: the surviving detections were not the objects + # the catalogue lists. + if len(objs) > 10: + rad0, _ = sep.flux_radius(sub, objs["x"], objs["y"], 6.0 * objs["a"], + 0.5, normflux=objs["flux"]) + fwhm0 = float(np.median(rad0) * 2.0) + min_npix = max(4, int(0.5 * np.pi * (fwhm0 / 2.0) ** 2)) + objs = objs[objs["npix"] >= min_npix] if len(objs) == 0: del data, sub, bkg return dict(sky=sky, rms=rms, fwhm=np.nan, nstars=0, From 2f268c7ca9e651118fbd2fde3c805cc9cc831f65 Mon Sep 17 00:00:00 2001 From: laurence Date: Tue, 21 Jul 2026 22:47:29 +0100 Subject: [PATCH 7/8] Add blind plate solving as the fallback when the seeded solver cannot match The seeded solver starts from the header's pointing, scale and roll angle and matches against Gaia. That works whenever the image and the catalogue hold recognisably the same stars, and fails when they do not: a single 60 second narrowband frame over a four degree field records a sparse, shallow population that overlaps poorly with any magnitude slice of Gaia, and no amount of window sliding fixes a population mismatch. A blind solver does not care about any of that. It builds geometric hashes from the image's own stars and looks them up in pre-built indexes, so it needs no pointing, no scale and no orientation - only pixels. blind.py drives nova.astrometry.net through astroquery, passing the plate scale as a hint where one is known, which turns a search over every possible scale into a search over one. Two properties worth stating plainly. It needs a free API key, and without one it says so clearly and the pipeline continues unsolved rather than failing. And it uploads the image to a third-party service, which is fine for these targets but is a fact worth knowing before pointing it at something unpublished. The result is verified rather than trusted. Whatever the service returns is matched back against Gaia locally and put through the same gate as the seeded solver - 40 stars and 1.5 px - so a solution has to be good, not merely returned. A wrong WCS remains worse than no WCS whoever produced it. --- pipeline/blind.py | 180 ++++++++++++++++++++++++++++++++++++++++++++++ pipeline/run.py | 11 +++ 2 files changed, 191 insertions(+) create mode 100644 pipeline/blind.py diff --git a/pipeline/blind.py b/pipeline/blind.py new file mode 100644 index 0000000..26812a2 --- /dev/null +++ b/pipeline/blind.py @@ -0,0 +1,180 @@ +"""Blind plate solving, for fields the seeded solver cannot handle. + +The seeded solver in `astrometry.py` starts from the header's pointing, scale +and roll angle and matches against Gaia. That works whenever the image and the +catalogue contain recognisably the same stars. It fails when they do not - a +single short narrowband frame over a four degree field records a sparse, +shallow star population that overlaps poorly with any magnitude slice of Gaia, +and no amount of window-sliding fixes a population mismatch. + +A blind solver does not care. It builds geometric hashes from the image's own +stars and looks them up in pre-built index files, so it needs no pointing, no +scale and no orientation - only the pixels. astrometry.net is the standard +implementation and this uses the hosted nova.astrometry.net service through +astroquery. + +Two things to know: + +* **It needs an API key.** Free, from an account at https://nova.astrometry.net + (Profile -> API Key). Provide it as the ASTROMETRY_API_KEY environment + variable, or in a file named `astrometry.key` beside the session or in + c:\\temp\\claudetemp. Without one this module reports that clearly and the + pipeline carries on unsolved rather than failing. +* **It uploads the image.** A stacked master goes to a third-party service. + That is fine for these targets and worth knowing before pointing it at + anything you would not publish. + +Results are still put through the same quality gate as the seeded solver: a +solution has to be good, not merely returned. +""" +import os + +import numpy as np +from astropy.io import fits +from astropy.wcs import WCS + +import layout + +KEY_FILES = ("astrometry.key", + r"c:\temp\claudetemp\astrometry.key") +MIN_STARS = 40 +MAX_RESID_PX = 1.5 + + +def find_key(session=None): + """The API key, from the environment or a key file.""" + key = os.environ.get("ASTROMETRY_API_KEY") + if key: + return key.strip() + candidates = list(KEY_FILES) + if session is not None: + candidates.insert(0, os.path.join(session.root, "astrometry.key")) + for path in candidates: + if os.path.exists(path): + with open(path, encoding="utf-8") as fh: + key = fh.read().strip() + if key: + return key + return None + + +def run(session, scale_hint=None, verbose=True): + """Solve the deepest master blind. Returns the same dict as the seeded + solver, or None.""" + key = find_key(session) + if not key: + if verbose: + print(" blind solve unavailable: no astrometry.net API key." + "\n Get one free at https://nova.astrometry.net " + "(Profile -> API Key), then set ASTROMETRY_API_KEY or write " + "it to c:\\temp\\claudetemp\\astrometry.key") + return None + + from astroquery.astrometry_net import AstrometryNet + + masters = os.path.join(session.root, layout.MASTERS) + deepest = session.filters[0] + path = os.path.join(masters, f"master-{deepest}.fit") + if not os.path.exists(path): + return None + + with fits.open(path) as hd: + image = hd[0].data.astype(np.float32) + header = hd[0].header + ny, nx = image.shape + + ast = AstrometryNet() + ast.api_key = key + ast.TIMEOUT = 600 + + # Give it the scale if it is known: a hint turns a search over every + # possible scale into a search over one, which is the difference between + # a minute and a quarter of an hour. + kwargs = dict(solve_timeout=600, publicly_visible="n") + if scale_hint: + kwargs.update(scale_units="arcsecperpix", scale_type="ev", + scale_est=float(scale_hint), scale_err=20) + if verbose: + print(f" blind solving {os.path.basename(path)} via " + f"nova.astrometry.net" + + (f" (scale hint {scale_hint:.3f} arcsec/px)" + if scale_hint else "")) + + try: + wcs_header = ast.solve_from_image(path, force_image_upload=True, + **kwargs) + except Exception as exc: # noqa: BLE001 + print(f" blind solve failed: {type(exc).__name__}: {exc}") + return None + if not wcs_header: + print(" blind solve returned no solution") + return None + + wcs = WCS(wcs_header) + cd = wcs.pixel_scale_matrix * 3600.0 + solved_scale = float(np.sqrt(abs(np.linalg.det(cd)))) + rot = float(np.degrees(np.arctan2(cd[0, 1], cd[1, 1]))) + centre = wcs.pixel_to_world(nx / 2.0, ny / 2.0) + + # Verify against Gaia rather than trusting the service, and put the answer + # through the same gate as the seeded solver. + import astrometry as seeded + from astropy import units as u + from astropy.coordinates import SkyCoord + from measure import measure_frame + + det = measure_frame(path, max_stars=1500) + src = det["xy"] + radius = 0.6 * np.hypot(nx, ny) * solved_scale / 3600.0 + cache = os.path.join(session.root, layout.INTERMEDIATES, "_gaia_blind.npz") + ra, dec, _ = seeded.fetch_catalogue(centre, radius, cache, limit=20000, + deep=True) + nmatch, resid = 0, None + if ra is not None and len(src) > 10: + world = SkyCoord(ra * u.deg, dec * u.deg) + px, py = wcs.world_to_pixel(world) + pairs = [] + for i, (cx, cy) in enumerate(zip(px, py)): + if not (0 <= cx < nx and 0 <= cy < ny): + continue + d = np.hypot(src[:, 0] - cx, src[:, 1] - cy) + j = int(np.argmin(d)) + if d[j] < 3.0: + pairs.append((i, j)) + if pairs: + ci = np.array([p[0] for p in pairs]) + ii = np.array([p[1] for p in pairs]) + qx, qy = wcs.world_to_pixel(world[ci]) + resid = np.hypot(qx - src[ii, 0], qy - src[ii, 1]) + nmatch = len(pairs) + + med = float(np.median(resid)) if resid is not None else float("nan") + if nmatch < MIN_STARS or not np.isfinite(med) or med > MAX_RESID_PX: + print(f" blind solution REJECTED on verification: {nmatch} " + f"Gaia stars at {med:.2f} px " + f"(need >= {MIN_STARS} and <= {MAX_RESID_PX})") + return None + + if verbose: + print(f" blind solve verified: {nmatch} Gaia stars, " + f"{med:.2f} px ({med * solved_scale:.2f} arcsec)") + print(f" centre {centre.to_string('hmsdms')}, " + f"{solved_scale:.4f} arcsec/px, PA {rot:.2f} deg") + + for filt in session.filters: + p = os.path.join(masters, f"master-{filt}.fit") + if not os.path.exists(p): + continue + with fits.open(p, mode="update") as hd: + for card in wcs.to_header().cards: + hd[0].header[card.keyword] = (card.value, card.comment) + hd[0].header["ASTRSOLV"] = ( + f"astrometry.net blind / {nmatch} Gaia stars / " + f"{med * solved_scale:.2f} arcsec", "blind plate solve") + hd.flush() + + return dict(nstars=nmatch, residual_px=med, + residual_arcsec=med * solved_scale, scale=solved_scale, + position_angle=rot, centre=centre.to_string("hmsdms"), + fov_arcmin=(nx * solved_scale / 60, ny * solved_scale / 60), + method="astrometry.net blind") diff --git a/pipeline/run.py b/pipeline/run.py index 3a985c5..ea9610b 100644 --- a/pipeline/run.py +++ b/pipeline/run.py @@ -84,6 +84,17 @@ def process(root, stages=STAGES, force=False, verbose=True): # An unsolved session still produces images; record and go on. print(f" solve failed: {type(exc).__name__}: {exc}") results["solve"] = None + if not results.get("solve"): + # The seeded solver needs the image and the catalogue to hold + # recognisably the same stars. When they do not, a blind solve + # ignores the pointing entirely and works from the pixels. + import blind + try: + results["solve"] = blind.run( + s, scale_hint=s.scale, verbose=verbose) + except Exception as exc: # noqa: BLE001 + print(f" blind solve failed: " + f"{type(exc).__name__}: {exc}") if "colour" in stages: import colour From 118dc36e926da1352578ab54442020371709e9e2 Mon Sep 17 00:00:00 2001 From: Laurence Date: Tue, 21 Jul 2026 23:00:23 +0100 Subject: [PATCH 8/8] solve: prefer the estate's own Astrometry.net over nova for blind solves The estate now runs Astrometry.net itself (ankh-morpork-infra astrometry/): solve-field against ~5 GB of local index files, sized from this repo's TELESCOPES.md so the whole iTelescope fleet's fields of view are covered. blind.py now tries it before nova.astrometry.net. Speed is the least of the reasons - a blind solve of a real DSS field returns in 0.7s where nova queues for minutes. The reasons that matter are that nova requires UPLOADING the master to a third party and holding an API key, and neither is necessary any more for the ordinary case. nova remains the fallback, so an estate outage costs speed and privacy rather than the ability to solve. pipeline/astrometry_net.py is the client and nothing more: stdlib-only POST, parses the returned .wcs so the full TAN solution is used rather than a re-derivation from the summary numbers. It deliberately does NOT decide whether to trust a solution - blind.py already verifies every blind solve against Gaia and applies a star-count and residual gate, and two gates that can disagree is worse than one that is trusted. The source ('estate' or 'nova') is threaded through the log lines, the returned method, and the ASTRSOLV card, because a year from now that card is the only way to tell whether a frame was solved in-house or uploaded. Also corrected astrometry.py's 'cannot solve without a blind solver' message, which has been untrue since run.py started falling through to blind.py. Tested against the live service: blind solve of a DSS2 field with a known centre returned within ~4 arcsec in 0.7s; the WCS round-trips through astropy; and an unreachable service falls through to nova instead of raising. --- pipeline/astrometry.py | 4 +- pipeline/astrometry_net.py | 142 +++++++++++++++++++++++++++++++++++++ pipeline/blind.py | 139 ++++++++++++++++++++++++++---------- state/DECISIONS.md | 23 ++++++ 4 files changed, 269 insertions(+), 39 deletions(-) create mode 100644 pipeline/astrometry_net.py diff --git a/pipeline/astrometry.py b/pipeline/astrometry.py index 4035d8a..5b29013 100644 --- a/pipeline/astrometry.py +++ b/pipeline/astrometry.py @@ -245,8 +245,10 @@ def run(session, verbose=True): centre = pointing(header) scale = plate_scale(header, session) if centre is None or not scale: + # Not fatal, and no longer the end of the road: returning None sends + # run.py to blind.py, which needs neither the pointing nor the scale. print(" no pointing or plate scale in the header - " - "cannot solve without a blind solver") + "handing off to the blind solver") return None # Match over the frame's INSCRIBED circle, not its circumscribed one. A # cone big enough to reach the corners also reaches well outside the diff --git a/pipeline/astrometry_net.py b/pipeline/astrometry_net.py new file mode 100644 index 0000000..1e14c0e --- /dev/null +++ b/pipeline/astrometry_net.py @@ -0,0 +1,142 @@ +"""Client for the estate's own Astrometry.net service. + + https://astrometry.ankh-morpork.discworld.network + +This is solve-field running on docker1/2/3 against ~5 GB of local index files. +It exists so that blind solving costs seconds instead of minutes and, more to +the point, so that no image has to leave the house to get solved. + +It is used by blind.py, which prefers this service and falls back to the hosted +nova.astrometry.net when it is unreachable. Everything a solution has to earn - +verification against Gaia, the star-count and residual gate, writing the WCS +into the masters - stays in blind.py, shared by both sources. This module only +gets a WCS out of the estate service; it does not decide whether to believe it. + +Only the stdlib is used for the request. A solve is one POST; pulling in +requests to do it would add a dependency for no gain. +""" +import base64 +import io +import os +import ssl +import time +import urllib.error +import urllib.request + +from astropy.io import fits +from astropy.wcs import WCS + +SERVICE = os.environ.get("ASTROMETRY_SERVICE", + "https://astrometry.ankh-morpork.discworld.network") +TIMEOUT = int(os.environ.get("ASTROMETRY_TIMEOUT", "180")) + + +def _context(): + """TLS context for the estate service. + + The service presents a step-ca certificate from the estate's internal CA. + On a machine that trusts the Discworld root CA this verifies normally. On + one that does not, verification fails - and the right fix is to install the + CA, not to turn checking off, so that is what the error message says. + ASTROMETRY_CA_BUNDLE points at the root CA if it is not in the system store. + ASTROMETRY_INSECURE=1 exists as an explicit, visible escape hatch; it is + never the default and it prints when used. + """ + bundle = os.environ.get("ASTROMETRY_CA_BUNDLE") + if os.environ.get("ASTROMETRY_INSECURE") == "1": + print(" WARNING: TLS verification disabled " + "(ASTROMETRY_INSECURE=1)") + ctx = ssl.create_default_context() + ctx.check_hostname = False + ctx.verify_mode = ssl.CERT_NONE + return ctx + return ssl.create_default_context(cafile=bundle) if bundle else \ + ssl.create_default_context() + + +def available(timeout=5): + """Is the service up and does it have indexes? Cheap enough to ask first.""" + try: + with urllib.request.urlopen(SERVICE + "/healthz", timeout=timeout, + context=_context()) as resp: + return resp.status == 200 + except Exception: # noqa: BLE001 + return False + + +def solve(path, centre=None, scale=None, radius_deg=2.0, timeout=TIMEOUT, + verbose=True): + """Solve one FITS file. Returns (WCS, info dict) or (None, reason). + + `centre` and `scale` are hints, not requirements - pass whatever the header + offered. A hinted solve turns a whole-sky search into a cone and is an order + of magnitude faster, but the whole point of this fallback is that it still + works when there is nothing to hint with. + """ + params = [] + if centre is not None: + params += [f"ra={centre.ra.deg:.6f}", f"dec={centre.dec.deg:.6f}", + f"radius={radius_deg:.3f}"] + if scale: + # +/-20% around the nominal scale: wide enough to absorb a focal + # reducer or a binning the header did not mention, tight enough to + # still rule out most index scales. + params += [f"scale_low={scale * 0.8:.4f}", + f"scale_high={scale * 1.2:.4f}"] + params.append(f"timeout={timeout}") + url = f"{SERVICE}/solve?" + "&".join(params) + + with open(path, "rb") as handle: + body = handle.read() + + started = time.time() + if verbose: + hint = "hinted" if (centre is not None or scale) else "blind" + print(f" astrometry.net ({hint}): " + f"{len(body) / 1e6:.1f} MB to {SERVICE}") + request = urllib.request.Request( + url, data=body, method="POST", + headers={"Content-Type": "application/octet-stream"}) + try: + with urllib.request.urlopen(request, timeout=timeout + 30, + context=_context()) as resp: + import json + payload = json.load(resp) + except urllib.error.HTTPError as exc: + # 422 is the service saying "no solution", which is an answer, not a + # fault. Anything else is the service or the network being broken, and + # the two should not read the same in the log. + import json + try: + detail = json.load(exc) + except Exception: # noqa: BLE001 + detail = {} + if exc.code == 422: + return None, "no solution found" + return None, f"service error HTTP {exc.code}: {detail.get('error', '')}" + except ssl.SSLCertVerificationError: + return None, ("TLS verification failed - install the Discworld root CA " + "or set ASTROMETRY_CA_BUNDLE to it") + except Exception as exc: # noqa: BLE001 + return None, f"{type(exc).__name__}: {exc}" + + if not payload.get("solved"): + return None, "no solution found" + + # The service hands back the .wcs file - a bare FITS header - base64 + # encoded. Parsing it here rather than trusting the summary numbers means + # the full distortion-free TAN solution is what gets written to the + # masters, not a re-derivation from ra/dec/scale. + with fits.open(io.BytesIO(base64.b64decode(payload["wcs_b64"]))) as hdul: + wcs = WCS(hdul[0].header) + + info = {k: payload.get(k) for k in + ("ra", "dec", "pixscale", "orientation", "parity", + "fieldw", "fieldh", "fieldunits")} + info["seconds"] = payload.get("seconds", round(time.time() - started, 1)) + if verbose: + print(f" astrometry.net solved in {info['seconds']}s: " + f"centre {info['ra']:.5f} {info['dec']:+.5f}, " + f"{info['pixscale']:.3f} arcsec/px, " + f"{info['fieldw']:.1f}' x {info['fieldh']:.1f}'") + return wcs, info diff --git a/pipeline/blind.py b/pipeline/blind.py index 26812a2..f397347 100644 --- a/pipeline/blind.py +++ b/pipeline/blind.py @@ -10,22 +10,33 @@ and no amount of window-sliding fixes a population mismatch. A blind solver does not care. It builds geometric hashes from the image's own stars and looks them up in pre-built index files, so it needs no pointing, no scale and no orientation - only the pixels. astrometry.net is the standard -implementation and this uses the hosted nova.astrometry.net service through -astroquery. +implementation, and there are now two of it available: -Two things to know: +1. **The estate's own** (`astrometry_net.py`) - solve-field on docker1/2/3 + against local index files. Tried first: the image never leaves the house, it + needs no API key, and it answers in seconds rather than queueing. Its index + set covers 5.6 arcmin to 33 degree quads, which is the whole iTelescope + fleet's range, but NOT arbitrarily narrow fields. +2. **The hosted nova.astrometry.net**, through astroquery. The fallback, for + when the estate service is unreachable or its index set cannot reach the + field. Two things to know about it: -* **It needs an API key.** Free, from an account at https://nova.astrometry.net - (Profile -> API Key). Provide it as the ASTROMETRY_API_KEY environment - variable, or in a file named `astrometry.key` beside the session or in - c:\\temp\\claudetemp. Without one this module reports that clearly and the - pipeline carries on unsolved rather than failing. -* **It uploads the image.** A stacked master goes to a third-party service. - That is fine for these targets and worth knowing before pointing it at - anything you would not publish. + * **It needs an API key.** Free, from an account at + https://nova.astrometry.net (Profile -> API Key). Provide it as the + ASTROMETRY_API_KEY environment variable, or in a file named + `astrometry.key` beside the session or in c:\\temp\\claudetemp. Without one + this module says so and the pipeline carries on unsolved rather than + failing. + * **It uploads the image.** A stacked master goes to a third-party service. + Fine for these targets, and worth knowing before pointing it at anything + you would not publish. This is the main reason the estate solver is tried + first. -Results are still put through the same quality gate as the seeded solver: a -solution has to be good, not merely returned. +Whichever source answers, the solution is put through the same quality gate as +the seeded solver - verified against Gaia, with a star count and residual it +has to meet. A solution has to be good, not merely returned. The source is +recorded in the ASTRSOLV card and in the returned `method`, so it is always +possible to tell afterwards whether a frame was solved in-house or uploaded. """ import os @@ -58,31 +69,50 @@ def find_key(session=None): return None -def run(session, scale_hint=None, verbose=True): - """Solve the deepest master blind. Returns the same dict as the seeded - solver, or None.""" +def _solve_via_estate(path, scale_hint=None, verbose=True): + """Try the estate's own solver first. Returns a WCS header, or None. + + Preferred over nova for three reasons, in order of how much they matter: + the image never leaves the house, it needs no API key, and it answers in + seconds rather than sitting in a public queue. If the service is down or + unreachable this returns None and the caller goes to nova, so the estate + being offline costs speed and privacy, not the ability to solve. + """ + try: + import astrometry_net + except ImportError: + return None + if not astrometry_net.available(): + if verbose: + print(f" estate solver unreachable " + f"({astrometry_net.SERVICE}) - falling back to nova") + return None + + wcs, info = astrometry_net.solve(path, scale=scale_hint, verbose=verbose) + if wcs is None: + # `info` is the reason string here. A genuine "no solution" is worth + # reporting but is not a reason to then try nova with the same pixels + # and the same index scales - nova's indexes are wider, so it is worth + # exactly one more attempt, and the caller makes it. + if verbose: + print(f" estate solver: {info}") + return None + return wcs.to_header() + + +def _solve_via_nova(path, session, scale_hint=None, verbose=True): + """The hosted service. Needs an API key and uploads the image.""" key = find_key(session) if not key: if verbose: - print(" blind solve unavailable: no astrometry.net API key." - "\n Get one free at https://nova.astrometry.net " - "(Profile -> API Key), then set ASTROMETRY_API_KEY or write " - "it to c:\\temp\\claudetemp\\astrometry.key") + print(" nova blind solve unavailable: no astrometry.net API" + " key.\n Get one free at https://nova.astrometry.net" + " (Profile -> API Key), then set ASTROMETRY_API_KEY or write" + " it to c:\\temp\\claudetemp\\astrometry.key") return None from astroquery.astrometry_net import AstrometryNet - masters = os.path.join(session.root, layout.MASTERS) - deepest = session.filters[0] - path = os.path.join(masters, f"master-{deepest}.fit") - if not os.path.exists(path): - return None - - with fits.open(path) as hd: - image = hd[0].data.astype(np.float32) - header = hd[0].header - ny, nx = image.shape - ast = AstrometryNet() ast.api_key = key ast.TIMEOUT = 600 @@ -104,10 +134,40 @@ def run(session, scale_hint=None, verbose=True): wcs_header = ast.solve_from_image(path, force_image_upload=True, **kwargs) except Exception as exc: # noqa: BLE001 - print(f" blind solve failed: {type(exc).__name__}: {exc}") + print(f" nova blind solve failed: {type(exc).__name__}: {exc}") return None if not wcs_header: - print(" blind solve returned no solution") + print(" nova blind solve returned no solution") + return None + return wcs_header + + +def run(session, scale_hint=None, verbose=True): + """Solve the deepest master blind. Returns the same dict as the seeded + solver, or None. + + Two possible sources, tried in that order: the estate's own solver, then + nova. Whichever answers, the solution then earns its place the same way - + verified against Gaia and put through the same gate as the seeded solver. + A blind solve that is merely *returned* is not a blind solve that is right. + """ + masters = os.path.join(session.root, layout.MASTERS) + deepest = session.filters[0] + path = os.path.join(masters, f"master-{deepest}.fit") + if not os.path.exists(path): + return None + + with fits.open(path) as hd: + image = hd[0].data.astype(np.float32) + ny, nx = image.shape + + source = "estate" + wcs_header = _solve_via_estate(path, scale_hint=scale_hint, verbose=verbose) + if wcs_header is None: + source = "nova" + wcs_header = _solve_via_nova(path, session, scale_hint=scale_hint, + verbose=verbose) + if wcs_header is None: return None wcs = WCS(wcs_header) @@ -150,13 +210,13 @@ def run(session, scale_hint=None, verbose=True): med = float(np.median(resid)) if resid is not None else float("nan") if nmatch < MIN_STARS or not np.isfinite(med) or med > MAX_RESID_PX: - print(f" blind solution REJECTED on verification: {nmatch} " - f"Gaia stars at {med:.2f} px " + print(f" {source} blind solution REJECTED on verification: " + f"{nmatch} Gaia stars at {med:.2f} px " f"(need >= {MIN_STARS} and <= {MAX_RESID_PX})") return None if verbose: - print(f" blind solve verified: {nmatch} Gaia stars, " + print(f" {source} blind solve verified: {nmatch} Gaia stars, " f"{med:.2f} px ({med * solved_scale:.2f} arcsec)") print(f" centre {centre.to_string('hmsdms')}, " f"{solved_scale:.4f} arcsec/px, PA {rot:.2f} deg") @@ -168,8 +228,11 @@ def run(session, scale_hint=None, verbose=True): with fits.open(p, mode="update") as hd: for card in wcs.to_header().cards: hd[0].header[card.keyword] = (card.value, card.comment) + # Name the SOURCE, not just "astrometry.net": a year from now the + # only way to tell whether a frame was solved in-house or uploaded + # to a public service is this card. hd[0].header["ASTRSOLV"] = ( - f"astrometry.net blind / {nmatch} Gaia stars / " + f"astrometry.net blind ({source}) / {nmatch} Gaia stars / " f"{med * solved_scale:.2f} arcsec", "blind plate solve") hd.flush() @@ -177,4 +240,4 @@ def run(session, scale_hint=None, verbose=True): residual_arcsec=med * solved_scale, scale=solved_scale, position_angle=rot, centre=centre.to_string("hmsdms"), fov_arcmin=(nx * solved_scale / 60, ny * solved_scale / 60), - method="astrometry.net blind") + method=f"astrometry.net blind ({source})") diff --git a/state/DECISIONS.md b/state/DECISIONS.md index cc15ca9..404eb1c 100644 --- a/state/DECISIONS.md +++ b/state/DECISIONS.md @@ -3,6 +3,29 @@ > A dated, append only log of decisions and their rationale. Newest at the top. Never > rewrite past entries; if a decision is reversed, add a new entry that says so. +## 2026-07-21: blind solving goes to the estate's own solver first, nova second + +The estate now runs Astrometry.net itself +(`https://astrometry.ankh-morpork.discworld.network`, config in +ankh-morpork-infra `astrometry/`): solve-field against ~5 GB of local index +files, sized from this repo's own TELESCOPES.md so it covers the whole +iTelescope fleet's fields of view (21.5 arcmin to 10.4 degrees). + +`blind.py` now tries it before nova.astrometry.net. The order is not about +speed, though it is far faster - a blind solve of a real DSS field came back in +0.7s against nova's minutes-in-a-queue. It is that **nova requires uploading +the image to a third party and an API key**, and neither is now necessary for +the ordinary case. nova stays as the fallback for when the estate service is +unreachable or the field is narrower than the local index set reaches, so +losing the estate costs speed and privacy, not the ability to solve at all. + +What deliberately did NOT change: the seeded solver in `astrometry.py` remains +the primary path (a good header beats a blind search), and every blind solution +- from either source - still has to pass the same Gaia verification and +star-count/residual gate. A wrong WCS is worse than no WCS, and a second solver +does not change that. The source is recorded in the ASTRSOLV card so it is +always possible to tell later which one produced a frame's astrometry. + ## 2026-07-21: the pipeline reads headers, not filenames The first pipeline parsed `calibrated-T32-...-Luminance-BIN2-W-300-001.fit` to