diff --git a/pipeline/astrometry.py b/pipeline/astrometry.py new file mode 100644 index 0000000..5b29013 --- /dev/null +++ b/pipeline/astrometry.py @@ -0,0 +1,424 @@ +"""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 +CATALOGUE_FETCH = 6000 # rich fields need far more than the match uses + + +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 _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) + 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: + # 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 - " + "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 + # 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") + + 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 * 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, + 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 + # 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. + 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: + 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)], max_control_points=60) + 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 match by pointing or asterism - solve failed") + return None + + 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(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 = [] + 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) + + # 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]))) + 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/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 new file mode 100644 index 0000000..f397347 --- /dev/null +++ b/pipeline/blind.py @@ -0,0 +1,243 @@ +"""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 there are now two of it available: + +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 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. + +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 + +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 _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(" 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 + + 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" nova blind solve failed: {type(exc).__name__}: {exc}") + return None + if not wcs_header: + 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) + 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" {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" {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") + + 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) + # 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 ({source}) / {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=f"astrometry.net blind ({source})") diff --git a/pipeline/colour.py b/pipeline/colour.py new file mode 100644 index 0000000..e2b1200 --- /dev/null +++ b/pipeline/colour.py @@ -0,0 +1,246 @@ +"""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 +from skimage.restoration import denoise_tv_chambolle + +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 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. + + 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 gradient: + for filt in list(data): + data[filt], _ = remove_gradient(data[filt], mask) + + if palette in PALETTE_MAP: + mapping = PALETTE_MAP[palette] + 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) + + refs = [float(np.percentile(c, 99.5)) for c in channels] + 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) + + 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: + lum = rgb.mean(axis=2) + + 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)") + + 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) + + 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 out + + +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/measure.py b/pipeline/measure.py new file mode 100644 index 0000000..a3b9b9e --- /dev/null +++ b/pipeline/measure.py @@ -0,0 +1,140 @@ +"""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=5, + deblend_cont=0.005) + # 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, + 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 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..ea9610b --- /dev/null +++ b/pipeline/run.py @@ -0,0 +1,181 @@ +"""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", + "science", "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 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 + 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") + 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()) 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")} 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 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).