Solve by pointing where the header allows it, and refuse doubtful solutions

Three fixes and one refusal, all found by working on the wide-field
narrowband session that would not solve.

The detection floor was a fixed 12 pixels, which quietly assumed a well
sampled star. At 0.53 arcsec/px with 5.9 px seeing that is right; at 3.5
arcsec/px the stars are undersampled at 1.75 px FWHM and cover a handful
of pixels each, so the floor discarded nearly every real star and kept
blends and galaxies instead. It now scales with the measured seeing, and
that session went from 305 usable detections to 600.

Matching now happens over the frame's inscribed circle rather than a cone
reaching its corners. On a 4 degree field the old cone covered 33 square
degrees of sky against 16 of image, so half the catalogue was not in the
picture at all.

Where the header records a roll angle - and iTelescope's does - the
orientation no longer has to be recovered from scratch. Rotation, scale
and parity are applied directly and only the residual pointing error is
searched, by histogramming every detection-to-catalogue offset and taking
the peak. Asterism matching remains as the fallback for headers that say
nothing about orientation.

The refusal matters most. With catalogue depth slices added, the wide
field produced a "solution" of 25 stars at 2.50 px residual, claiming a
centre half a degree from the pointing. A genuine solve on that field
matches hundreds of stars and refines to well under a pixel. Accepting it
would have silently corrupted every position in the science catalogue, so
a solution must now reach 40 stars AND 1.5 px or it is discarded and the
session is reported as unsolved. A wrong WCS is worse than no WCS.

NGC 2030 improves from 0.31 to 0.12 arcsec on 123 stars under the
pointing-based match. NGC 2070 is honestly unsolved.

pipeline/astrometry.py

@@ -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,6 +310,9 @@ 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.
+    if best == "pointing":
+        pass
+    else:
       best = None
       for offset in (0, 15, 40, 80):
         for n in (120, 200, 60):
@@ -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,9 +340,10 @@ 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
 
+    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)))
@@ -232,7 +355,7 @@ def run(session, verbose=True):
     # 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])))

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
+    # Saturated cores and cosmic-ray hits make poor registration anchors and
-    # and poor seeing estimates, so they are excluded before anything is
+    # poor seeing estimates, so they go first.
-    # measured from them.
+    objs = objs[(objs["flag"] == 0) & (objs["npix"] < 3000)]
-    objs = objs[(objs["flag"] == 0) & (objs["npix"] > 12) &
+
-                (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,