Fix the blend filter's FOV pre-filter units and bound its memory

RMS/Astrometry/StarFilters.py

@@ -13,15 +13,67 @@
 
 import numpy as np
 
+from scipy.spatial import cKDTree
+
 from RMS.Astrometry.StarClasses import PairedStars
-from RMS.Astrometry.ApplyAstrometry import extinctionCorrectionTrueToApparent, raDecToXYPP
+from RMS.Astrometry.ApplyAstrometry import extinctionCorrectionTrueToApparent, raDecToXYPP, \
+    getFOVSelectionRadius, xyToRaDecPP
 
 
 # Default filtering parameters
 DEFAULT_PHOTOMETRIC_SIGMA = 2.5
 DEFAULT_BLEND_FWHM_MULT = 2.0  # Multiplier of FWHM for blending detection radius
 DEFAULT_BLEND_MAG_MARGIN = 0.3  # Margin above limiting magnitude for blend check
 
+# Margin applied to the FOV selection radius when pre-filtering the catalog. The radius
+#   already circumscribes the image corners, so the margin only has to cover neighbours
+#   sitting a few pixels outside the frame - it is generous rather than tight.
+DEFAULT_FOV_RADIUS_MARGIN = 1.5
+
+
+def catalogStarsInFOV(catalog_ra, catalog_dec, platepar, jd, margin=DEFAULT_FOV_RADIUS_MARGIN):
+    """ Mask of the catalog stars which lie inside the FOV cone at the given time.
+
+    Stars behind the camera can reverse-project into valid-looking pixel coordinates, so a
+    cone around the pointing direction is used to reject them before projecting.
+
+    Arguments:
+        catalog_ra: [ndarray] Catalog star right ascensions (deg).
+        catalog_dec: [ndarray] Catalog star declinations (deg).
+        platepar: [Platepar] Platepar for the FOV geometry.
+        jd: [float] Julian date.
+
+    Keyword arguments:
+        margin: [float] Multiplier applied to the FOV selection radius.
+            Default is DEFAULT_FOV_RADIUS_MARGIN.
+
+    Returns:
+        in_fov: [ndarray] Boolean mask, True for stars inside the cone.
+    """
+    # Radius which includes the image corners, computed by projecting them through the
+    #   platepar (distortion included) instead of assuming the central F_scale holds all
+    #   the way out to the corners
+    fov_radius = min(getFOVSelectionRadius(platepar)*margin, 90)
+
+    # Centre the cone on the pointing at THIS jd. platepar.RA_d/dec_d is the pointing at
+    #   platepar.JD, and on an alt-az camera the two drift ~15 deg/hour apart, which would
+    #   eat the margin above whenever a platepar is reused across a night.
+    _, ra_centre, dec_centre, _ = xyToRaDecPP([jd], [platepar.X_res/2.0], [platepar.Y_res/2.0], [1],
+                                              platepar, extinction_correction=False, jd_time=True)
+
+    ra_centre = np.radians(ra_centre[0])
+    dec_centre = np.radians(dec_centre[0])
+
+    ra_rad = np.radians(catalog_ra)
+    dec_rad = np.radians(catalog_dec)
+
+    # Spherical angular distance from the pointing direction to each catalog star
+    cos_ang_dist = (np.sin(dec_centre)*np.sin(dec_rad)
+                    + np.cos(dec_centre)*np.cos(dec_rad)*np.cos(ra_rad - ra_centre))
+    cos_ang_dist = np.clip(cos_ang_dist, -1, 1)
+
+    return np.degrees(np.arccos(cos_ang_dist)) < fov_radius
+
 
 def filterPhotometricOutliers(paired_stars, platepar, jd, sigma_threshold=DEFAULT_PHOTOMETRIC_SIGMA,
                                verbose=False):
@@ -150,23 +202,7 @@ def filterBlendedStars(paired_stars, catalog_stars, platepar, jd, lim_mag,
     # Filter to stars actually in front of the camera (within FOV + margin)
     # This prevents false positives from stars behind the camera that could
     # project to valid-looking pixel coordinates
-    ra_rad = np.radians(catalog_ra)
-    dec_rad = np.radians(catalog_dec)
-    ra_center = np.radians(platepar.RA_d)
-    dec_center = np.radians(platepar.dec_d)
-
-    # Spherical angular distance from camera pointing to each catalog star
-    cos_ang_dist = (np.sin(dec_center) * np.sin(dec_rad) +
-                    np.cos(dec_center) * np.cos(dec_rad) * np.cos(ra_rad - ra_center))
-    cos_ang_dist = np.clip(cos_ang_dist, -1, 1)
-    ang_dist_deg = np.degrees(np.arccos(cos_ang_dist))
-
-    # Estimate FOV radius from platepar (diagonal / 2 * scale, with margin)
-    fov_diagonal = np.sqrt(platepar.X_res**2 + platepar.Y_res**2)
-    fov_radius = (fov_diagonal / 2) * platepar.F_scale * 1.5  # 50% margin
-    fov_radius = min(fov_radius, 90)  # Cap at 90 degrees
-
-    in_fov = ang_dist_deg < fov_radius
+    in_fov = catalogStarsInFOV(catalog_ra, catalog_dec, platepar, jd)
     catalog_ra = catalog_ra[in_fov]
     catalog_dec = catalog_dec[in_fov]
 
@@ -198,15 +234,35 @@ def filterBlendedStars(paired_stars, catalog_stars, platepar, jd, lim_mag,
             np.array(matched_ra_list), np.array(matched_dec_list), jd, platepar)
         blend_radii = np.array(blend_radii)
 
-        # Compute distance from each matched star to all bright catalog stars using broadcasting
-        # Shape: (n_matched, n_catalog)
-        dx = all_matched_x[:, np.newaxis] - catalog_x[np.newaxis, :]
-        dy = all_matched_y[:, np.newaxis] - catalog_y[np.newaxis, :]
-        dist_matrix = np.sqrt(dx**2 + dy**2)
-
-        # Check for neighbors within each star's blend radius (excluding self)
-        has_neighbor = np.any(
-            (dist_matrix < blend_radii[:, np.newaxis]) & (dist_matrix > 0.1), axis=1)
+        # Every matched star is inside the image, so a catalog star further outside the
+        # frame than the largest blend radius can never be a blend neighbour. This is
+        # exact in the units that matter (pixels) and removes most of what the cone
+        # pre-filter necessarily lets through.
+        edge_margin = np.max(blend_radii)
+        in_img = ((catalog_x > -edge_margin) & (catalog_x < platepar.X_res + edge_margin)
+                  & (catalog_y > -edge_margin) & (catalog_y < platepar.Y_res + edge_margin))
+        catalog_x = catalog_x[in_img]
+        catalog_y = catalog_y[in_img]
+
+        if len(catalog_x) == 0:
+            return paired_stars, 0
+
+        # Nearest-neighbour lookup rather than an (n_matched x n_catalog) distance matrix.
+        # O(N log M) with constant memory, so a deep catalog cannot allocate the multi-GB
+        # matrices that used to get the process OOM-killed, and there is no chunk size to
+        # tune. Same reasoning as RMS/Astrometry/MatchStars.py.
+        tree = cKDTree(np.column_stack([catalog_x, catalog_y]))
+
+        # A matched star normally coincides with its own catalog entry (d ~ 0), which the
+        # d > 0.1 test below drops. Ask for a few neighbours rather than one so the search
+        # still reaches a genuine neighbour when the catalog holds duplicate entries at
+        # that position. Missing neighbours come back as inf, which fails both tests.
+        n_neighbors = min(3, len(catalog_x))
+        nn_dist, _ = tree.query(np.column_stack([all_matched_x, all_matched_y]), k=n_neighbors)
+        nn_dist = np.reshape(nn_dist, (len(check_indices), n_neighbors))
+
+        # Neighbours within each star's own blend radius (excluding self)
+        has_neighbor = np.any((nn_dist < blend_radii[:, np.newaxis]) & (nn_dist > 0.1), axis=1)
 
         for k, idx in enumerate(check_indices):
             if has_neighbor[k]: