diff --git a/pipeline/astrometry.py b/pipeline/astrometry.py deleted file mode 100644 index 5b29013..0000000 --- a/pipeline/astrometry.py +++ /dev/null @@ -1,424 +0,0 @@ -"""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 deleted file mode 100644 index 1e14c0e..0000000 --- a/pipeline/astrometry_net.py +++ /dev/null @@ -1,142 +0,0 @@ -"""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 deleted file mode 100644 index f397347..0000000 --- a/pipeline/blind.py +++ /dev/null @@ -1,243 +0,0 @@ -"""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 deleted file mode 100644 index e2b1200..0000000 --- a/pipeline/colour.py +++ /dev/null @@ -1,246 +0,0 @@ -"""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 deleted file mode 100644 index cdde6f8..0000000 --- a/pipeline/finals.py +++ /dev/null @@ -1,257 +0,0 @@ -"""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 deleted file mode 100644 index a3b9b9e..0000000 --- a/pipeline/measure.py +++ /dev/null @@ -1,140 +0,0 @@ -"""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 deleted file mode 100644 index c4a522d..0000000 --- a/pipeline/register.py +++ /dev/null @@ -1,159 +0,0 @@ -"""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 deleted file mode 100644 index 2b1b025..0000000 --- a/pipeline/report.py +++ /dev/null @@ -1,188 +0,0 @@ -"""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 deleted file mode 100644 index ea9610b..0000000 --- a/pipeline/run.py +++ /dev/null @@ -1,181 +0,0 @@ -"""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 deleted file mode 100644 index d6c4e82..0000000 --- a/pipeline/science.py +++ /dev/null @@ -1,297 +0,0 @@ -"""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 404eb1c..cc15ca9 100644 --- a/state/DECISIONS.md +++ b/state/DECISIONS.md @@ -3,29 +3,6 @@ > 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 8342390..efbcbf9 100644 --- a/state/TODO.md +++ b/state/TODO.md @@ -18,40 +18,8 @@ one telescope. Branch `generic-pipeline`. filter grouping, palette detection (LRGB / RGB / SHO / HOO / MONO), manifest writing. Tested against all five archived sessions. -**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. +**Next:** `measure` -> `register/stack` -> `solve` -> `compose` -> `report`, +then a `run.py` orchestrator with resumable stages. **The blocker before any of that:** the disk is full (see below).