diff --git a/pipeline/astrometry.py b/pipeline/astrometry.py index 4887980..4035d8a 100644 --- a/pipeline/astrometry.py +++ b/pipeline/astrometry.py @@ -35,6 +35,7 @@ import layout import measure as measure_mod SEED_STARS = 600 +CATALOGUE_FETCH = 6000 # rich fields need far more than the match uses def pointing(header): @@ -140,6 +141,94 @@ def project(ra, dec, centre, scale, parity): dy.to_value(u.arcsec) / scale]) + +def _roll_angle(header): + """Camera angle from the header, if the telescope recorded it.""" + for key in ("HIERARCH RollAngle", "RollAngle", "POSANG", "CROTA2"): + val = header.get(key) + if val is None: + continue + try: + return float(val) + except (TypeError, ValueError): + continue + return None + + +def match_by_pointing(src, ra, dec, centre, scale, roll, shape, verbose=True): + """Pair detections with catalogue stars using a known camera angle. + + Asterism matching exists to recover an UNKNOWN orientation. When the header + records the roll angle - and iTelescope's does - the whole geometry is + already known bar a small pointing error, and a direct search is both more + robust and far faster. Rotation, scale and parity are applied, then the + residual translation is found by histogramming every detection-to-catalogue + offset and taking the peak: a real solution puts thousands of pairs in one + bin, and noise spreads flat. + """ + ny, nx = shape + # Slide down the catalogue's magnitude ranking as well as trying the four + # geometric conventions. In a field as rich as the LMC the brightest few + # hundred catalogue stars are far brighter than anything a 60 second + # narrowband frame records, so the two "brightest N" lists describe + # different populations and overlap barely at all. + best = None + slices = [(0, 400), (200, 800), (600, 1400), (1200, 2400), (2000, 4000)] + for lo, hi in slices: + if lo >= len(ra): + break + ra_s, dec_s = ra[lo:hi], dec[lo:hi] + for parity in (1.0, -1.0): + for sign in (1.0, -1.0): + th = np.radians(sign * roll) + cat = project(ra_s, dec_s, centre, scale, parity) + rot = np.column_stack([ + cat[:, 0] * np.cos(th) - cat[:, 1] * np.sin(th), + cat[:, 0] * np.sin(th) + cat[:, 1] * np.cos(th)]) + cat_px = rot + np.array([nx / 2.0, ny / 2.0]) + + # Every pairwise offset, histogrammed. Limited to the brightest of + # each list to keep this to a few million comparisons. + a = src[:400] + b = cat_px[:400] + dx = (a[:, None, 0] - b[None, :, 0]).ravel() + dy = (a[:, None, 1] - b[None, :, 1]).ravel() + lim = 0.25 * max(nx, ny) + keep = (np.abs(dx) < lim) & (np.abs(dy) < lim) + if keep.sum() < 50: + continue + bins = np.arange(-lim, lim + 12, 12) + H, xe, ye = np.histogram2d(dx[keep], dy[keep], bins=[bins, bins]) + iy, ix = np.unravel_index(np.argmax(H), H.shape) + peak = H[iy, ix] + if best is None or peak > best[0]: + off = (0.5 * (xe[iy] + xe[iy + 1]), + 0.5 * (ye[ix] + ye[ix + 1])) + best = (peak, parity, sign, off, cat_px, lo, hi) + if best is None: + return None + peak, parity, sign, off, cat_px, lo, hi = best + if verbose: + print(f" pointing match: catalogue[{lo}:{hi}], parity " + f"{parity:+.0f}, roll {sign:+.0f}, offset " + f"({off[0]:+.0f}, {off[1]:+.0f}) px, peak {int(peak)} pairs") + if peak < 12: + return None + shifted = cat_px + np.array(off) + idx_base = lo + pairs = [] + for i, (cx, cy) in enumerate(shifted): + if not (0 <= cx < nx and 0 <= cy < ny): + continue + d = np.hypot(src[:, 0] - cx, src[:, 1] - cy) + j = int(np.argmin(d)) + if d[j] < 18: + pairs.append((i + idx_base, j)) + if len(pairs) < 10: + return None + return np.array([p[0] for p in pairs]), np.array([p[1] for p in pairs]) + + def run(session, verbose=True): """Solve the deepest master; write the WCS into every master.""" masters_dir = os.path.join(session.root, layout.MASTERS) @@ -159,7 +248,14 @@ def run(session, verbose=True): print(" no pointing or plate scale in the header - " "cannot solve without a blind solver") return None - radius = 1.15 * 0.5 * np.hypot(nx, ny) * scale / 3600.0 + # Match over the frame's INSCRIBED circle, not its circumscribed one. A + # cone big enough to reach the corners also reaches well outside the + # frame - on a 4 degree field that is 33 square degrees of catalogue + # against 16 of image, so half the catalogue's brightest stars are not in + # the picture at all and the two "brightest N" lists barely overlap. + # The inscribed circle is inside the frame whatever the camera angle, so + # both lists then describe the same patch of sky. + radius = 0.95 * 0.5 * min(nx, ny) * scale / 3600.0 if verbose: print(f" pointing {centre.to_string('hmsdms')}, " f"{scale:.4f} arcsec/px, search {radius:.3f} deg") @@ -167,21 +263,44 @@ def run(session, verbose=True): xy, _ = measure_mod.load(session) # Detect on the master itself: it is deeper than any single frame. from measure import measure_frame - det = measure_frame(path, max_stars=SEED_STARS) + det = measure_frame(path, max_stars=SEED_STARS * 3) src = det["xy"] + # Same restriction on the detections, so the two lists cover one region. + r_pix = 0.95 * 0.5 * min(nx, ny) + inside = np.hypot(src[:, 0] - nx / 2.0, src[:, 1] - ny / 2.0) < r_pix + src = src[inside][:SEED_STARS] if len(src) < 12: print(f" only {len(src)} stars in the master - too few") return None cache = os.path.join(session.root, layout.INTERMEDIATES, "_gaia.npz") os.makedirs(os.path.dirname(cache), exist_ok=True) - ra, dec, gmag = fetch_catalogue(centre, radius, cache) + ra, dec, gmag = fetch_catalogue(centre, radius, cache, + limit=CATALOGUE_FETCH) if ra is None: print(" no catalogue available - skipping the solve") return None if verbose: print(f" {len(ra)} catalogue stars, {len(src)} detected") + # If the header recorded the camera angle, use it: a direct geometric + # match is far more reliable than recovering the orientation from scratch. + roll = _roll_angle(header) + seed = None + if roll is not None: + if verbose: + print(f" header roll angle {roll:.2f} deg") + seed = match_by_pointing(src, ra, dec, centre, scale, roll, + (ny, nx), verbose=verbose) + if seed is not None: + ci, ii = seed + world0 = SkyCoord(ra[ci] * u.deg, dec[ci] * u.deg) + wcs = fit_wcs_from_points((src[ii, 0], src[ii, 1]), world0, + proj_point="center", projection="TAN") + best = "pointing" + else: + best = None + # Matching needs the two lists to contain the SAME stars, and "brightest" # does not guarantee that. The detector deliberately rejects saturated # cores, so on a well exposed frame the brightest detections are not the @@ -191,8 +310,11 @@ def run(session, verbose=True): # # So slide a window down the catalogue's magnitude ranking as well as # varying the sample size, and take the best match found. - best = None - for offset in (0, 15, 40, 80): + if best == "pointing": + pass + else: + best = None + for offset in (0, 15, 40, 80): for n in (120, 200, 60): hi = min(offset + n, len(ra)) if hi - offset < 25 or len(src) < 25: @@ -203,7 +325,7 @@ def run(session, verbose=True): parity) try: tform, (cat_m, img_m) = aa.find_transform( - cat_xy, src[:max(n, 120)]) + cat_xy, src[:max(n, 120)], max_control_points=60) except Exception: # noqa: BLE001 continue if best is None or len(cat_m) > len(best[0]): @@ -218,21 +340,22 @@ def run(session, verbose=True): break if best is None: - print(" no asterism match in either parity - solve failed") + print(" no match by pointing or asterism - solve failed") return None - cat_m, img_m, parity, _offset = best - cat_xy = project(ra, dec, centre, scale, parity) - idx = [int(np.argmin(np.hypot(cat_xy[:, 0] - px, cat_xy[:, 1] - py))) - for px, py in cat_m] - world = SkyCoord(ra[idx] * u.deg, dec[idx] * u.deg) - wcs = fit_wcs_from_points((img_m[:, 0], img_m[:, 1]), world, - proj_point="center", projection="TAN") + if best != "pointing": + cat_m, img_m, parity, _offset = best + cat_xy = project(ra, dec, centre, scale, parity) + idx = [int(np.argmin(np.hypot(cat_xy[:, 0] - px, cat_xy[:, 1] - py))) + for px, py in cat_m] + world = SkyCoord(ra[idx] * u.deg, dec[idx] * u.deg) + wcs = fit_wcs_from_points((img_m[:, 0], img_m[:, 1]), world, + proj_point="center", projection="TAN") # Refine: the seed match is a handful of stars, so use the approximate # solution to pair every catalogue source with its nearest detection. all_world = SkyCoord(ra * u.deg, dec * u.deg) - nmatch, resid = len(cat_m), None + nmatch, resid = (len(seed[0]) if seed is not None else len(cat_m)), None for tol in (4.0, 2.0): px, py = wcs.world_to_pixel(all_world) pairs = [] @@ -253,6 +376,19 @@ def run(session, verbose=True): resid = np.hypot(qx - src[ii, 0], qy - src[ii, 1]) nmatch = len(pairs) + # A quality gate, because a WRONG solution is worse than none: every + # position in the science catalogue would inherit the error, silently. A + # genuine solve on a field like this matches hundreds of stars and refines + # to well under a pixel; a spurious alignment matches a couple of dozen and + # stalls at several. Both conditions must hold. + MIN_STARS, MAX_RESID_PX = 40, 1.5 + if resid is None or nmatch < MIN_STARS or np.median(resid) > MAX_RESID_PX: + got = "n/a" if resid is None else f"{np.median(resid):.2f} px" + print(f" REJECTED: {nmatch} stars at {got} residual " + f"(need >= {MIN_STARS} stars and <= {MAX_RESID_PX} px). " + f"Treating as unsolved rather than trusting a doubtful WCS.") + return None + cd = wcs.pixel_scale_matrix * 3600.0 solved_scale = float(np.sqrt(abs(np.linalg.det(cd)))) rot = float(np.degrees(np.arctan2(cd[0, 1], cd[1, 1]))) diff --git a/pipeline/measure.py b/pipeline/measure.py index 2593ee1..a3b9b9e 100644 --- a/pipeline/measure.py +++ b/pipeline/measure.py @@ -32,13 +32,26 @@ def measure_frame(path, max_stars=MAX_STARS): rms = float(bkg.globalrms) sub = data - bkg.back() - objs = sep.extract(sub, DETECT_SIGMA, err=rms, minarea=9, + objs = sep.extract(sub, DETECT_SIGMA, err=rms, minarea=5, deblend_cont=0.005) - # Saturated cores and cosmic-ray hits both make poor registration anchors - # and poor seeing estimates, so they are excluded before anything is - # measured from them. - objs = objs[(objs["flag"] == 0) & (objs["npix"] > 12) & - (objs["npix"] < 3000)] + # Saturated cores and cosmic-ray hits make poor registration anchors and + # poor seeing estimates, so they go first. + objs = objs[(objs["flag"] == 0) & (objs["npix"] < 3000)] + + # The minimum size has to follow the SAMPLING, not be a fixed number of + # pixels. A fixed floor of 12 px assumes a well sampled star: at 0.53 + # arcsec/px with 5.9 px seeing that is right, but at 3.5 arcsec/px the + # stars are undersampled at 1.75 px FWHM and cover only a handful of + # pixels each, so the same floor discards nearly every real star and keeps + # blends and galaxies instead. That is what made the wide-field sessions + # impossible to plate solve: the surviving detections were not the objects + # the catalogue lists. + if len(objs) > 10: + rad0, _ = sep.flux_radius(sub, objs["x"], objs["y"], 6.0 * objs["a"], + 0.5, normflux=objs["flux"]) + fwhm0 = float(np.median(rad0) * 2.0) + min_npix = max(4, int(0.5 * np.pi * (fwhm0 / 2.0) ** 2)) + objs = objs[objs["npix"] >= min_npix] if len(objs) == 0: del data, sub, bkg return dict(sky=sky, rms=rms, fwhm=np.nan, nstars=0,