From cd62e52284f54714ce23e41a1b4c0617b72014ab Mon Sep 17 00:00:00 2001 From: Luc Busquin <133058544+Cybis320@users.noreply.github.com> Date: Fri, 24 Jul 2026 09:56:49 -0700 Subject: [PATCH 1/5] Fix the blend filter's FOV pre-filter units and bound its memory The FOV radius pre-filter multiplied the pixel diagonal by F_scale (px/deg) instead of dividing, yielding a huge value that always hit the 90-degree cap - the filter was a silent no-op, so a deep catalog fed millions of stars into the all-pairs blend matrices and the process got OOM-killed. Same units bug fixed in SkyFit2's copy. The pairwise distance computation now also runs in bounded-memory catalog chunks, so peak memory stays flat regardless of catalog depth. No behavior change: the same stars are flagged as blended. --- RMS/Astrometry/StarFilters.py | 31 ++++++++++++++++++++----------- Utils/SkyFit2.py | 3 ++- 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/RMS/Astrometry/StarFilters.py b/RMS/Astrometry/StarFilters.py index 0f4daedfe..eee75fceb 100644 --- a/RMS/Astrometry/StarFilters.py +++ b/RMS/Astrometry/StarFilters.py @@ -161,9 +161,9 @@ def filterBlendedStars(paired_stars, catalog_stars, platepar, jd, lim_mag, 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) + # Estimate FOV radius from platepar (F_scale is px/deg), 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 = (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 @@ -198,15 +198,24 @@ 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) + # Compute distance from each matched star to all bright catalog stars using + # broadcasting, in catalog chunks so peak memory stays bounded no matter how + # many catalog stars survived the pre-filters (a deep catalog fed through the + # broken FOV pre-filter above used to allocate multi-GB matrices here and get + # the process OOM-killed) + # Shape per chunk: (n_matched, chunk) + n_matched = len(check_indices) + chunk_size = max(1, int(5e6) // n_matched) + has_neighbor = np.zeros(n_matched, dtype=bool) + for c0 in range(0, len(catalog_x), chunk_size): + c1 = c0 + chunk_size + dx = all_matched_x[:, np.newaxis] - catalog_x[np.newaxis, c0:c1] + dy = all_matched_y[:, np.newaxis] - catalog_y[np.newaxis, c0:c1] + 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) for k, idx in enumerate(check_indices): if has_neighbor[k]: diff --git a/Utils/SkyFit2.py b/Utils/SkyFit2.py index a42764d5a..5110feed0 100644 --- a/Utils/SkyFit2.py +++ b/Utils/SkyFit2.py @@ -5895,8 +5895,9 @@ def count_matches_at_lm(test_lm): ang_dist_deg = np.degrees(np.arccos(cos_ang_dist)) # FOV radius with margin (stars behind camera have ang_dist > 90) + # F_scale is px/deg, so divide to convert the pixel diagonal to degrees fov_diagonal = np.sqrt(self.platepar.X_res**2 + self.platepar.Y_res**2) - fov_radius = (fov_diagonal / 2) * self.platepar.F_scale * 1.5 + fov_radius = (fov_diagonal / 2) / self.platepar.F_scale * 1.5 fov_radius = min(fov_radius, 90) in_fov = ang_dist_deg < fov_radius From 1910164c90e40885b260739d26d87c2dc478748b Mon Sep 17 00:00:00 2001 From: Luc Busquin <133058544+Cybis320@users.noreply.github.com> Date: Tue, 28 Jul 2026 09:30:12 -0700 Subject: [PATCH 2/5] Use the exact FOV selection radius and centre the catalog cone on the image time --- RMS/Astrometry/StarFilters.py | 71 +++++++++++++++++++++++++---------- Utils/SkyFit2.py | 24 ++---------- 2 files changed, 55 insertions(+), 40 deletions(-) diff --git a/RMS/Astrometry/StarFilters.py b/RMS/Astrometry/StarFilters.py index eee75fceb..5bffc0273 100644 --- a/RMS/Astrometry/StarFilters.py +++ b/RMS/Astrometry/StarFilters.py @@ -12,9 +12,9 @@ from __future__ import print_function, division, absolute_import import numpy as np - 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 @@ -22,6 +22,55 @@ 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 +199,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 (F_scale is px/deg), 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] diff --git a/Utils/SkyFit2.py b/Utils/SkyFit2.py index 5110feed0..0ef7b7840 100644 --- a/Utils/SkyFit2.py +++ b/Utils/SkyFit2.py @@ -154,7 +154,8 @@ def computeSolarSystemMagnitude(body_name, body, sun, time): limitingMagnitude, screenNudgeToAzAltDelta, fovCentreZenithDirection from RMS.Astrometry.AtmosphericExtinction import atmosphericExtinctionCorrection from RMS.Astrometry.StarClasses import CatalogStar, GeoPoint, PlanetPoint, PairedStars -from RMS.Astrometry.StarFilters import filterPhotometricOutliers, filterBlendedStars +from RMS.Astrometry.StarFilters import filterPhotometricOutliers, filterBlendedStars, \ + catalogStarsInFOV from RMS.Astrometry.Conversions import date2JD, JD2HourAngle, trueRaDec2ApparentAltAz, \ apparentAltAz2TrueRADec, J2000_JD, jd2Date, datetime2JD, JD2LST, geo2Cartesian, vector2RaDec, raDec2Vector from RMS.Astrometry.AstrometryNet import astrometryNetSolve @@ -5881,26 +5882,7 @@ def count_matches_at_lm(test_lm): # Filter catalog to stars actually in front of camera (prevent back-projection) # Use angular distance in celestial coordinates, not self.filterCatalogStarsInsideFOV # which incorrectly uses self.cat_lim_mag instead of the test catalog's LM - catalog_ra = test_catalog[:, 0] - catalog_dec = test_catalog[:, 1] - ra_rad = np.radians(catalog_ra) - dec_rad = np.radians(catalog_dec) - ra_center = np.radians(self.platepar.RA_d) - dec_center = np.radians(self.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)) - - # FOV radius with margin (stars behind camera have ang_dist > 90) - # F_scale is px/deg, so divide to convert the pixel diagonal to degrees - fov_diagonal = np.sqrt(self.platepar.X_res**2 + self.platepar.Y_res**2) - fov_radius = (fov_diagonal / 2) / self.platepar.F_scale * 1.5 - fov_radius = min(fov_radius, 90) - - in_fov = ang_dist_deg < fov_radius + in_fov = catalogStarsInFOV(test_catalog[:, 0], test_catalog[:, 1], self.platepar, jd) test_catalog = test_catalog[in_fov] if len(test_catalog) == 0: From e41c4adf4daf6eb8c71894d0b91dba1f0b67f9c9 Mon Sep 17 00:00:00 2001 From: Luc Busquin <133058544+Cybis320@users.noreply.github.com> Date: Tue, 28 Jul 2026 09:30:12 -0700 Subject: [PATCH 3/5] Replace the blend distance matrix with an image-bounds cull and a KD-tree --- RMS/Astrometry/StarFilters.py | 50 ++++++++++++++++++++++------------- 1 file changed, 32 insertions(+), 18 deletions(-) diff --git a/RMS/Astrometry/StarFilters.py b/RMS/Astrometry/StarFilters.py index 5bffc0273..40d68fcb0 100644 --- a/RMS/Astrometry/StarFilters.py +++ b/RMS/Astrometry/StarFilters.py @@ -12,6 +12,9 @@ from __future__ import print_function, division, absolute_import import numpy as np + +from scipy.spatial import cKDTree + from RMS.Astrometry.StarClasses import PairedStars from RMS.Astrometry.ApplyAstrometry import extinctionCorrectionTrueToApparent, raDecToXYPP, \ getFOVSelectionRadius, xyToRaDecPP @@ -231,24 +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, in catalog chunks so peak memory stays bounded no matter how - # many catalog stars survived the pre-filters (a deep catalog fed through the - # broken FOV pre-filter above used to allocate multi-GB matrices here and get - # the process OOM-killed) - # Shape per chunk: (n_matched, chunk) - n_matched = len(check_indices) - chunk_size = max(1, int(5e6) // n_matched) - has_neighbor = np.zeros(n_matched, dtype=bool) - for c0 in range(0, len(catalog_x), chunk_size): - c1 = c0 + chunk_size - dx = all_matched_x[:, np.newaxis] - catalog_x[np.newaxis, c0:c1] - dy = all_matched_y[:, np.newaxis] - catalog_y[np.newaxis, c0:c1] - 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]: From b812ff94237fd891df69f8448751a644e606ee0f Mon Sep 17 00:00:00 2001 From: Luc Busquin <133058544+Cybis320@users.noreply.github.com> Date: Tue, 28 Jul 2026 09:30:12 -0700 Subject: [PATCH 4/5] Add tests for the blended-star catalog pre-filters --- Tests/TestStarFilters.py | 317 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 317 insertions(+) create mode 100644 Tests/TestStarFilters.py diff --git a/Tests/TestStarFilters.py b/Tests/TestStarFilters.py new file mode 100644 index 000000000..07751da1a --- /dev/null +++ b/Tests/TestStarFilters.py @@ -0,0 +1,317 @@ +""" Tests for the catalog pre-filters used by the blended-star rejection. + +Covers the FOV cone pre-filter (units, exactness and the epoch of its centre) and the +neighbour search in filterBlendedStars. +""" + +from __future__ import print_function, division, absolute_import + +import os + +import pytest + +np = pytest.importorskip("numpy") + +from RMS.Astrometry.ApplyAstrometry import getFOVSelectionRadius, raDecToXYPP, xyToRaDecPP +from RMS.Astrometry.StarClasses import CatalogStar, PairedStars +from RMS.Astrometry.StarFilters import catalogStarsInFOV, filterBlendedStars +from RMS.Formats.Platepar import Platepar + + +# Real 720p platepar shipped with RMS. Using a fitted platepar rather than a hand-made one +# matters here: it carries actual distortion coefficients, and its reference RA/Dec is +# consistent with the projection, which a synthetic platepar is not. +TEMPLATE = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), + 'share', 'platepar_templates', 'template_generic_720p_4mm.cal') + + +def makePlatepar(): + """ Platepar from the shipped 4 mm 720p template. + + The template's reference JD is left alone - its RA/Dec, Alt/Az and pointing are only + mutually consistent at the epoch it was fitted for. + """ + + if not os.path.isfile(TEMPLATE): + pytest.skip("platepar template not available: {:s}".format(TEMPLATE)) + + pp = Platepar() + pp.read(TEMPLATE) + pp.refraction = False + + return pp + + +def pixelToRaDec(pp, x, y, jd=None): + """ RA/Dec of the given image coordinates, as arrays. """ + + if jd is None: + jd = pp.JD + + x = np.atleast_1d(np.asarray(x, dtype=np.float64)) + y = np.atleast_1d(np.asarray(y, dtype=np.float64)) + + _, ra, dec, _ = xyToRaDecPP(len(x)*[jd], x, y, np.ones(len(x)), pp, + extinction_correction=False, jd_time=True) + + return ra, dec + + +def fovCentre(pp, jd=None): + """ RA/Dec of the FOV centre at the given time. """ + + ra, dec = pixelToRaDec(pp, pp.X_res/2.0, pp.Y_res/2.0, jd=jd) + + return ra[0], dec[0] + + +def angularSeparationDeg(ra1, dec1, ra2, dec2): + """ Great-circle separation between two points (deg). """ + + ra1, dec1, ra2, dec2 = map(np.radians, (ra1, dec1, ra2, dec2)) + cos_sep = np.sin(dec1)*np.sin(dec2) + np.cos(dec1)*np.cos(dec2)*np.cos(ra1 - ra2) + + return np.degrees(np.arccos(np.clip(cos_sep, -1, 1))) + + +class TestCatalogStarsInFOV(object): + + def test_units_regression(self): + """ The cone must not be a hemisphere pass-through. + + Guards the original bug: multiplying by F_scale (px/deg) instead of dividing put + the radius in the thousands of degrees, so the min(..., 90) cap always fired and + every catalog star in front of the camera survived. + """ + + pp = makePlatepar() + ra_c, dec_c = fovCentre(pp) + + # A star 85 deg off-axis is in front of the camera but nowhere near the FOV, so it + # is only accepted if the cone has degenerated to a hemisphere + far_ra = np.array([ra_c]) + far_dec = np.array([dec_c - 85.0]) + + assert angularSeparationDeg(ra_c, dec_c, far_ra[0], far_dec[0]) > 80 + assert not catalogStarsInFOV(far_ra, far_dec, pp, pp.JD)[0] + + def test_image_corners_are_inside(self): + """ Every pixel of the image, corners included, must survive the pre-filter. """ + + pp = makePlatepar() + + xs = [0, pp.X_res, 0, pp.X_res, pp.X_res/2.0] + ys = [0, pp.Y_res, pp.Y_res, 0, pp.Y_res/2.0] + ra, dec = pixelToRaDec(pp, xs, ys) + + assert np.all(catalogStarsInFOV(ra, dec, pp, pp.JD)) + + # Still true with essentially no margin - getFOVSelectionRadius circumscribes the + # image, so the corners land on the radius itself. A hair over 1.0 keeps this off + # an exact float boundary, since the two sides compute the separation differently. + assert np.all(catalogStarsInFOV(ra, dec, pp, pp.JD, margin=1.001)) + + def test_radius_is_at_least_the_selection_radius(self): + """ The accepted cone must cover getFOVSelectionRadius, i.e. the whole image. + + Direct regression guard on the units bug, phrased as the property that matters: + anything the platepar can image has to be considered. + """ + + pp = makePlatepar() + ra_c, dec_c = fovCentre(pp) + fov_radius = getFOVSelectionRadius(pp) + + # Walk outwards along a meridian and find where the filter starts rejecting + offsets = np.arange(0, 90, 0.5) + ra = np.full(offsets.shape, ra_c) + dec = dec_c - offsets + + accepted = catalogStarsInFOV(ra, dec, pp, pp.JD) + seps = angularSeparationDeg(ra_c, dec_c, ra, dec) + + assert np.max(seps[accepted]) >= fov_radius + + def test_cone_centre_follows_jd(self): + """ The cone is centred on the pointing at jd, not on platepar.RA_d/dec_d. + + An alt-az camera sweeps ~15 deg/hour in RA, so a platepar reused across a night + must not keep filtering against the pointing it was fitted at. + """ + + pp = makePlatepar() + + jd_later = pp.JD + 0.5 # 12 h later, i.e. roughly the opposite side of the sky + ra_later, dec_later = fovCentre(pp, jd=jd_later) + ra_arr, dec_arr = np.array([ra_later]), np.array([dec_later]) + + # Where the camera is actually pointing at jd_later + assert catalogStarsInFOV(ra_arr, dec_arr, pp, jd_later)[0] + + # ...is not where it pointed at platepar.JD + assert not catalogStarsInFOV(ra_arr, dec_arr, pp, pp.JD)[0] + + def test_empty_catalog(self): + """ An empty catalog gives an empty mask rather than an error. """ + + pp = makePlatepar() + mask = catalogStarsInFOV(np.array([]), np.array([]), pp, pp.JD) + + assert len(mask) == 0 + + +def buildCatalog(pp, pixel_positions, mag=5.0, jd=None): + """ Catalog array [ra, dec, mag] for stars at the given image coordinates. """ + + xs = [p[0] for p in pixel_positions] + ys = [p[1] for p in pixel_positions] + ra, dec = pixelToRaDec(pp, xs, ys, jd=jd) + + return np.column_stack([ra, dec, np.full(len(ra), mag)]) + + +def buildPairedStars(pp, pixel_positions, fwhm=3.0, mag=5.0, jd=None): + """ PairedStars whose catalog coordinates project back to the given image coordinates. """ + + catalog = buildCatalog(pp, pixel_positions, mag=mag, jd=jd) + + paired = PairedStars() + for (x, y), (ra, dec, star_mag) in zip(pixel_positions, catalog): + paired.addPair(x, y, fwhm, 1000.0, CatalogStar(ra, dec, star_mag)) + + return paired + + +class TestFilterBlendedStars(object): + + # 5 matched stars spread across the frame - the filter no-ops below 5 pairs + BASE = [(300, 200), (500, 250), (700, 300), (900, 350), (400, 500)] + + def test_close_neighbour_is_flagged(self): + """ A catalog star at 1.5x FWHM is inside the 2x FWHM blend radius. """ + + pp = makePlatepar() + paired = buildPairedStars(pp, self.BASE, fwhm=3.0) + + # 4.5 px away from the first matched star, blend radius is 2*3 = 6 px + catalog = buildCatalog(pp, self.BASE + [(304.5, 200)]) + + filtered, removed = filterBlendedStars(paired, catalog, pp, pp.JD, 6.0) + + assert removed == 1 + assert len(filtered) == len(paired) - 1 + + def test_distant_neighbour_is_kept(self): + """ A catalog star at 3x FWHM is outside the blend radius. """ + + pp = makePlatepar() + paired = buildPairedStars(pp, self.BASE, fwhm=3.0) + + # 9 px away, blend radius is 6 px + catalog = buildCatalog(pp, self.BASE + [(309, 200)]) + + filtered, removed = filterBlendedStars(paired, catalog, pp, pp.JD, 6.0) + + assert removed == 0 + assert len(filtered) == len(paired) + + def test_self_match_is_not_a_blend(self): + """ A star's own catalog entry sits at d ~ 0 and must not flag it. """ + + pp = makePlatepar() + paired = buildPairedStars(pp, self.BASE, fwhm=3.0) + catalog = buildCatalog(pp, self.BASE) + + _, removed = filterBlendedStars(paired, catalog, pp, pp.JD, 6.0) + + assert removed == 0 + + def test_duplicate_catalog_entry_does_not_hide_a_blend(self): + """ Duplicate entries at the star's own position must not mask a real neighbour. + + The nearest-neighbour search asks for several neighbours precisely so that two + coincident catalog rows cannot crowd out the genuine one. + """ + + pp = makePlatepar() + paired = buildPairedStars(pp, self.BASE, fwhm=3.0) + + # The first matched star's position twice, plus a genuine neighbour at 4.5 px + catalog = buildCatalog(pp, self.BASE + [(300, 200), (304.5, 200)]) + + _, removed = filterBlendedStars(paired, catalog, pp, pp.JD, 6.0) + + assert removed == 1 + + def test_matches_bruteforce_reference(self): + """ The KD-tree search must agree with an explicit O(N*M) distance scan. + + Replaces the chunked-vs-unchunked check the chunking approach needed. + """ + + pp = makePlatepar() + rng = np.random.RandomState(42) + + matched_positions = [(float(x), float(y)) for x, y in + rng.uniform([100, 100], [1100, 600], size=(40, 2))] + catalog_positions = [(float(x), float(y)) for x, y in + rng.uniform([100, 100], [1100, 600], size=(400, 2))] + + fwhm = 4.0 + paired = buildPairedStars(pp, matched_positions, fwhm=fwhm) + catalog = buildCatalog(pp, matched_positions + catalog_positions) + + filtered, removed = filterBlendedStars(paired, catalog, pp, pp.JD, 6.0) + + # Reference: project everything the same way the filter does, then scan + cat_x, cat_y = raDecToXYPP(catalog[:, 0], catalog[:, 1], pp.JD, pp) + matched_x = np.array([p[0] for p in matched_positions]) + matched_y = np.array([p[1] for p in matched_positions]) + + blend_radius = 2.0*fwhm + expected = 0 + for mx, my in zip(matched_x, matched_y): + dist = np.hypot(cat_x - mx, cat_y - my) + if np.any((dist < blend_radius) & (dist > 0.1)): + expected += 1 + + assert expected > 0, "test data should contain at least one blend" + assert removed == expected + assert len(filtered) == len(matched_positions) - expected + + def test_deep_catalog_outside_the_fov_is_cheap(self): + """ A catalog covering the whole sky must not be carried into the distance search. + + The pre-filters exist so a deep catalog cannot allocate an (n_matched x n_catalog) + matrix; this checks the result is unaffected by stars that cannot possibly blend. + """ + + pp = makePlatepar() + paired = buildPairedStars(pp, self.BASE, fwhm=3.0) + + near_catalog = buildCatalog(pp, self.BASE + [(304.5, 200)]) + + # Whole-sky catalog of equally bright stars, almost all of it behind the camera or + # far off-axis + rng = np.random.RandomState(7) + ra_all = rng.uniform(0, 360, 20000) + dec_all = np.degrees(np.arcsin(rng.uniform(-1, 1, 20000))) + sky = np.column_stack([ra_all, dec_all, np.full(len(ra_all), 5.0)]) + + catalog = np.vstack([near_catalog, sky]) + + _, removed = filterBlendedStars(paired, catalog, pp, pp.JD, 6.0) + + assert removed == 1 + + def test_below_minimum_pairs_is_a_noop(self): + """ Fewer than 5 pairs returns the input untouched. """ + + pp = makePlatepar() + paired = buildPairedStars(pp, self.BASE[:4], fwhm=3.0) + catalog = buildCatalog(pp, self.BASE[:4]) + + filtered, removed = filterBlendedStars(paired, catalog, pp, pp.JD, 6.0) + + assert removed == 0 + assert filtered is paired From 608a5c898694fa2939b64b98498e2f8effcd9170 Mon Sep 17 00:00:00 2001 From: Luc Busquin <133058544+Cybis320@users.noreply.github.com> Date: Tue, 28 Jul 2026 10:07:18 -0700 Subject: [PATCH 5/5] Reuse filterCatalogStarsInsideFOV for the LM search catalog cone --- Utils/SkyFit2.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/Utils/SkyFit2.py b/Utils/SkyFit2.py index 0ef7b7840..8b0372bae 100644 --- a/Utils/SkyFit2.py +++ b/Utils/SkyFit2.py @@ -154,8 +154,7 @@ def computeSolarSystemMagnitude(body_name, body, sun, time): limitingMagnitude, screenNudgeToAzAltDelta, fovCentreZenithDirection from RMS.Astrometry.AtmosphericExtinction import atmosphericExtinctionCorrection from RMS.Astrometry.StarClasses import CatalogStar, GeoPoint, PlanetPoint, PairedStars -from RMS.Astrometry.StarFilters import filterPhotometricOutliers, filterBlendedStars, \ - catalogStarsInFOV +from RMS.Astrometry.StarFilters import filterPhotometricOutliers, filterBlendedStars from RMS.Astrometry.Conversions import date2JD, JD2HourAngle, trueRaDec2ApparentAltAz, \ apparentAltAz2TrueRADec, J2000_JD, jd2Date, datetime2JD, JD2LST, geo2Cartesian, vector2RaDec, raDec2Vector from RMS.Astrometry.AstrometryNet import astrometryNetSolve @@ -5879,13 +5878,17 @@ def count_matches_at_lm(test_lm): if test_catalog is None or len(test_catalog) == 0: return 0, 0 - # Filter catalog to stars actually in front of camera (prevent back-projection) - # Use angular distance in celestial coordinates, not self.filterCatalogStarsInsideFOV - # which incorrectly uses self.cat_lim_mag instead of the test catalog's LM - in_fov = catalogStarsInFOV(test_catalog[:, 0], test_catalog[:, 1], self.platepar, jd) - test_catalog = test_catalog[in_fov] - - if len(test_catalog) == 0: + # Filter the catalog to stars actually in front of the camera, which prevents the + # back-projection fold. filterCatalogStarsInsideFOV takes an explicit lim_mag, so + # the test catalog's own LM is applied rather than self.cat_lim_mag, and it centres + # the cone via computeCentreRADec() at the CURRENT image time. Centring on + # platepar.RA_d/dec_d instead is only valid at platepar.JD: after a whole-night + # refit that epoch can sit hours away from the displayed frame, and the cone then + # points at the wrong sky, so the true FOV stars are excluded and only + # back-projection folds survive. + _, test_catalog = self.filterCatalogStarsInsideFOV(test_catalog, lim_mag=test_lm) + + if (test_catalog is None) or (len(test_catalog) == 0): return 0, 0 cat_x, cat_y, _ = getCatalogStarsImagePositions(test_catalog, jd, self.platepar)