Skip to content

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

Open
Cybis320 wants to merge 5 commits into
prereleasefrom
fix-blend-filter-fov-units
Open

Fix the blend filter's FOV pre-filter units and bound its memory#939
Cybis320 wants to merge 5 commits into
prereleasefrom
fix-blend-filter-fov-units

Conversation

@Cybis320

@Cybis320 Cybis320 commented Jul 24, 2026

Copy link
Copy Markdown
Contributor

The FOV radius pre-filter in filterBlendedStars multiplied the pixel diagonal by F_scale (px/deg) instead of dividing, producing a huge value that always hit the 90-degree cap — the pre-filter was a silent no-op. With a deep catalog, every catalog star then flowed into the all-pairs blend matrices, allocating multi-GB arrays and getting the process OOM-killed. The same units bug was fixed in SkyFit2's copy of the FOV check.

Both copies are now a single helper, catalogStarsInFOV() in RMS/Astrometry/StarFilters.py, which uses getFOVSelectionRadius() — the established call elsewhere in the codebase, exact because it projects the image corners through the platepar with distortion included — and centres the cone on the pointing at the image's own jd rather than at platepar.JD. On an alt-az camera those drift ~15 deg/hour apart, so a platepar reused across a night was filtering against a stale pointing.

The distance search is now a cKDTree nearest-neighbour lookup instead of an (n_matched × n_catalog) matrix: O(N log M), constant memory, no chunk size to tune. Catalog stars are also culled to the image bounds plus the largest blend radius before the search — every matched star is inside the frame, so anything further out than that cannot blend with one.

Behavior change: the units fix is not neutral. It narrows the cone from a 90-degree hemisphere to the actual FOV plus margin, so catalog stars far off-axis are no longer considered — precisely the stars the pre-filter's comment says it exists to exclude, since reverse projection can fold them back into valid-looking pixel coordinates. That is a false-positive reduction, and blended-star counts may differ. The chunking-to-KD-tree change is behaviour-preserving.

Known related issue, not fixed here: the magnitude pre-filter (StarFilters.py:140-141) is a silent no-op from the same root cause — both callers pass the catalog's own limiting magnitude as lim_mag, so bright_mask selects mag < catalog_limit + 0.3, i.e. everything. Fixing it changes which stars are rejected as blended, so it wants its own PR and its own before/after on station data. Noted here so a future deep-catalog OOM isn't read as "the FOV fix didn't work".

Note on reading the diff: much of it is re-indentation rather than semantic change.

Tests: Tests/TestStarFilters.py, 12 cases covering the cone (units, corner coverage, epoch of the centre) and the neighbour search (blend radius boundaries, self-match, duplicate catalog entries, agreement with an explicit O(N×M) reference scan).

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.

@dvida dvida left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Verdict: the units fix is correct and should land

Confirmed against the tree: platepar.F_scale is documented px/deg (RMS/Formats/Platepar.py:229, and the .cal loader converts arcsec/px → px/deg at Platepar.py:2223-2225). Multiplying was dimensionally wrong, so min(fov_radius, 90) always fired and the pre-filter was a hemisphere pass-through — exactly as described.

The fix also errs in the safe direction, which is worth stating since the whole point of the pre-filter is not to drop real stars: F_scale is the central scale, and for both gnomonic (narrow-field, where arctan compresses the edges) and equidistant-fisheye lenses the linear estimate (diag/2)/F_scale is ≥ the true corner angle. So no in-FOV catalog star is newly excluded.

The chunking is functionally correct too — has_neighbor |= ORs across chunks and the trailing partial chunk is handled by the slice.

Findings below. Inline comments cover four; the rest are on lines outside the diff hunks, so they are here. None of them block the units fix.


1. The magnitude pre-filter is also a silent no-op — same root cause, not fixed

RMS/Astrometry/StarFilters.py:140-141:

max_mag = lim_mag + mag_margin
bright_mask = catalog_stars[:, 2] < max_mag

Both call sites pass the catalog's own limiting magnitude as lim_mag:

  • RMS/Astrometry/AutoPlatepar.py:657-658 passes config.catalog_mag_limit, and the catalog_stars it forwards is the whole-sky catalog from loadCatalogStars(config, config.catalog_mag_limit) (AutoPlatepar.py:784) — never FOV-subset.
  • Utils/SkyFit2.py:8620-8628 passes self.cat_lim_mag, which is the catalog load depth (SkyFit2.py:2483), and the LM search pushes it deep at SkyFit2.py:5312.

So bright_mask selects mag < catalog_limit + 0.3 — i.e. everything. The docstring at StarFilters.py:123 says lim_mag is the "Current limiting magnitude for star detection", which is not what either caller supplies.

This is the co-equal half of why the matrices blew up: two pre-filters, both no-ops, one fixed here. Fixing it is out of scope for this PR, but worth calling out in the description so the next deep-catalog OOM doesn't get read as "the FOV fix didn't work".

2. FOV centre comes from platepar.RA_d/dec_d while the distances are for an arbitrary jd

StarFilters.py:155-156 and SkyFit2.py (same block). RA_d is the apparent RA of the image centre at platepar.JD. For an alt-az fixed camera the centre's RA drifts ~15°/hour, so a platepar reused across a night puts the cone centre tens of degrees off the true pointing — which eats most of the margin this PR just introduced.

SkyFit2.py:5832, a few dozen lines above the block you patched, already does it the robust way:

_, ra_c, dec_c, _ = xyToRaDecPP([jd2Date(jd)], [pp.X_res/2], [pp.Y_res/2], [1], pp,
                                extinction_correction=False)

In AutoPlatepar this is currently benign — the platepar is fitted at the same jd (AutoPlatepar.py:601) right before the filter call (:657) — but it's load-bearing on an assumption nothing enforces.

3. Missing image-bounds cull after projection — one line, bigger win than the cone filter

Right after raDecToXYPP at StarFilters.py:177: every matched star is inside the image, and the blend radius is a few px, so a catalog star projecting outside the frame can never be a blend neighbour. Culling

in_img = (catalog_x > -m) & (catalog_x < platepar.X_res + m) \
       & (catalog_y > -m) & (catalog_y < platepar.Y_res + m)

with m ≈ a few px cuts the working set roughly another order of magnitude and makes the chunk-size question mostly academic.

4. "No behavior change" isn't quite right

The chunking half is behaviour-preserving. The units half is not — it narrows the cone from a 90° hemisphere to ~76° on a typical 1280×720 / ~15 px/deg setup, so catalog stars 76–90° off-axis are no longer considered. Those are precisely the stars the comment at StarFilters.py:150-152 says the check exists to exclude: far off-axis stars the reverse projection can fold back into valid-looking pixel coordinates.

That's a false-positive reduction — a real and desirable change. Better to say that than "No behavior change", so nobody is surprised when the blended-star count differs.

5. No test committed

The description says "verified chunked == unchunked on test data" but nothing ships. Tests/ already has the pattern (Tests/TestMatchStars.py, Tests/TestMath.py). A Tests/TestStarFilters.py with three cases would lock all of this in:

  • fov_radiusgetFOVSelectionRadius(platepar) for a synthetic platepar — fails on the multiply version, passes here. Direct regression guard on the units bug.
  • chunked == unchunked with the chunk size forced to 1, 3, n-1, n, n+1 against a catalog length that is not a multiple of the chunk size.
  • a star with a neighbour at 1.5×FWHM is flagged; one at 3×FWHM is not.

Minimum I'd want before merge: keep the units fix, correct the description per #4, and tighten the chunk budget (see the inline comment on chunk_size). The getFOVSelectionRadius and cKDTree points are both net line deletions, so they're cheap if you want them in the same PR.

Comment thread RMS/Astrometry/StarFilters.py Outdated
# 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

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The units are now right, but this hand-rolled estimate doesn't need to exist — getFOVSelectionRadius(platepar) (RMS/Astrometry/ApplyAstrometry.py:500-531) projects the four image corners through the actual platepar, distortion included, and returns the max angular separation from centre. StarFilters.py:17 already imports from that module.

It's the established call everywhere else in the codebase, including in SkyFit2 itself at lines 4559, 5826 and 12325, plus CheckFit.py:133, ApplyRecalibrate.py:224, NNalign.py:133 and AddCelestialGrid.py. NNalign.py:136 even uses the exact idiom you'd want here:

fov_radius = getFOVSelectionRadius(platepar)
fov_radius_margin = fov_radius * 1.5

Three reasons this is more than style:

  1. It's exact, so the margin can shrink. On a 1280×720 / ~15 px/deg camera this formula gives ≈76° where the true corner radius is ≈50°. Cone solid angle ∝ (1−cos r): 0.76 vs 0.36. Swapping in the helper roughly halves the surviving catalog again, on top of what this PR already recovers — which is the actual goal here.
  2. It handles lenses the linear approximation only approximates. F_scale is the central scale; the corner angle depends on the distortion polynomial, which getFOVSelectionRadius evaluates rather than assumes.
  3. It removes the duplication that caused this bug. This ~20-line cone block is byte-identical to Utils/SkyFit2.py:5899-5921, and this PR patches both copies by hand — the same failure mode that let the units error sit in two places. Extracting one shared helper (or just calling getFOVSelectionRadius from both) makes the third copy impossible.

Comment thread RMS/Astrometry/StarFilters.py Outdated
# the process OOM-killed)
# Shape per chunk: (n_matched, chunk)
n_matched = len(check_indices)
chunk_size = max(1, int(5e6) // n_matched)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The budget counts elements, not bytes, and doesn't account for how many arrays are live at once. 5e6 float64 = 40 MB per array, but np.sqrt(dx**2 + dy**2) on the next lines keeps ~6 full-size temporaries alive simultaneously — dx, dy, dx**2, dy**2, their sum, and the sqrt result — so real peak transient is ~240 MB, plus the two bool masks. RMS runs on 2 GB Raspberry Pi 4s, so "bounded" is bounded at a figure that can still hurt.

Two cheap tightenings:

  • Drop the budget to ~5e5 elements (~30 MB peak). Loop overhead is noise next to the arithmetic.
  • Compare squared distances and skip sqrt entirely: d2 = dx*dx + dy*dy, then (d2 < blend_radii[:, None]**2) & (d2 > 0.01). One fewer full-size temporary and no transcendental.

Nit while you're here: int(5e6) inside a floor-division reads like a unit slip waiting to happen. Write it as a named module-level constant next to the existing DEFAULT_* values at StarFilters.py:21-23, with the byte budget in the comment rather than the element count.

Comment thread RMS/Astrometry/StarFilters.py Outdated
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)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Chunking bounds the memory but leaves the time at O(n_matched × n_catalog) — still ~10⁸ distance evaluations per call with a deep catalog even after the FOV fix, and filterBlendedStars runs inside the recalibration loop.

The codebase already has the answer and documents the reasoning: RMS/Astrometry/MatchStars.py:2 opens with "scipy.spatial.cKDTree for O(N log M) performance instead of O(NM)"*. cKDTree is imported in Utils/SkyFit2.py:127, RMS/ExtractStars.py:31 and RMS/Formats/Platepar.py:40, and scipy>=1.0.0 is in requirements.txt:13.

tree = cKDTree(np.column_stack([catalog_x, catalog_y]))
# k=2: the star itself (d~0) plus its nearest neighbour
d, _ = tree.query(np.column_stack([all_matched_x, all_matched_y]), k=2)
has_neighbor = (d[:, 1] < blend_radii) & (d[:, 1] > 0.1)

Constant memory, no chunk-size knob to tune, and it deletes the loop rather than sizing it. k=2 plus > 0.1 reproduces the existing self-exclusion; if a matched star can legitimately have two catalog entries at d≈0, use query_ball_point(..., r=blend_radii.max()) and filter per-star instead.

Comment thread Utils/SkyFit2.py Outdated
# 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

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same fix, correct — but this is the duplication half of the point I left on StarFilters.py:166. This block is byte-identical to RMS/Astrometry/StarFilters.py:150-171, and both copies had to be patched by hand in this PR. Meanwhile this same file already calls getFOVSelectionRadius(self.platepar) at lines 4559, 5826 and 12325 — including at 5826, inside _computeSeasonalStarVariation, ~70 lines above this one.

Swapping both copies for getFOVSelectionRadius(platepar) * 1.5 makes a third divergent copy impossible and is a net line deletion.

@Cybis320

Copy link
Copy Markdown
Contributor Author

Thanks — took findings 2, 3 and 5 plus both inline suggestions, in 1910164c..b812ff94. One correction on the arithmetic behind the getFOVSelectionRadius swap, because it points the opposite way from the review and it changes why the swap is worth doing.

The "halves the catalog again" number doesn't hold

The comparison is 76° with the 1.5 margin against ~50° without it, so it's measuring the margin, not the exactness. Keeping the * 1.5 your own snippet has, the swap is not a reduction at all — and on a wide lens it goes the other way, because the linear estimate under-estimates the corner radius rather than over-estimating it.

Measured on the two platepars shipped in share/platepar_templates/:

template getFOVSelectionRadius (diag/2)/F_scale linear error
generic_720p_4mm (15.6 px/deg) 55.73° 47.03° −15.6%
generic_720p_6mm (24.1 px/deg) 30.95° 30.47° −1.6%

So on the 4 mm lens the exact cone is wider: 83.6° vs 70.6° at margin 1.5. Against a uniform whole-sky catalog, as a fraction of the hemisphere the original let through — 1.00 for the old capped-at-90 behaviour, 0.67 for this PR's linear form, 0.89 with the helper, 0.44 for the helper with no margin at all.

That also means the margin is doing real work, not just padding: at margin 1.0 the linear estimate (47.03°) is smaller than the true corner radius (55.73°), so it would clip catalog stars that are genuinely in frame. Tests/TestStarFilters.py::test_image_corners_are_inside fails against the linear formula for exactly that reason.

The helper still belongs here — it's exact by construction, it removes the duplicated block, and it's what the rest of the codebase calls. It just isn't a size win, so I'd rather it land on those grounds than on a number that won't show up in the logs.

2 — cone centre now follows jd (1910164c)

Taken as suggested, using the xyToRaDecPP idiom from SkyFit2.py:5832. This is the one finding that could clip real stars, so it seemed worth doing here rather than deferring.

Both copies of the cone are now one helper, catalogStarsInFOV(catalog_ra, catalog_dec, platepar, jd, margin=1.5) in StarFilters.py, called from filterBlendedStars and from SkyFit2.count_matches_at_lm. Net deletion, and a third divergent copy is no longer possible.

3 — image-bounds cull (e41c4adf)

Added, with m = the largest blend radius rather than a constant, since that's exactly the distance past the frame edge at which a catalog star stops being able to blend with anything.

Measured smaller than you estimated, though: 4.2×, not an order of magnitude. Reverse projection folds a lot of far off-axis sky back into valid pixel coordinates — 24% of what survives the cone still lands within the frame bounds. Which is the same effect that motivates the cone in the first place, so the two guards are genuinely complementary; neither replaces the other.

Chunk budget → replaced rather than tuned (e41c4adf)

Went with cKDTree, so the element budget, the int(5e6), and the chunk loop are all gone rather than resized. O(N log M), constant memory, nothing to tune.

One deviation from your snippet: k=3 instead of k=2, with both/all columns tested rather than just d[:, 1]. With k=2 and a fixed d[:, 1], two coincident catalog rows at a matched star's own position crowd out the real neighbour and the blend is missed — and if the star's own catalog entry isn't in the filtered subset, d[:, 1] skips past a genuine neighbour sitting at d[:, 0]. Testing every returned distance against > 0.1 handles both. test_duplicate_catalog_entry_does_not_hide_a_blend covers the first case.

5 — tests (b812ff94)

Tests/TestStarFilters.py, 12 cases, following the Tests/TestMatchStars.py pattern. They're built on share/platepar_templates/template_generic_720p_4mm.cal rather than a hand-made platepar — a synthetic one has no distortion coefficients (so the helper and the linear estimate agree exactly, and the whole question becomes invisible) and its RA_d/dec_d isn't consistent with its own projection, which quietly breaks any test about the cone centre.

Run against the three versions of the pre-filter:

                                                units_regression  corners_inside  centre_follows_jd
original (F_scale multiply, platepar centre)          FAIL             pass             FAIL
this PR at cd62e522 (linear, platepar centre)         pass             FAIL             FAIL
with 1910164c..b812ff94                               pass             pass             pass

The blend-search half is checked against an explicit O(N×M) reference scan on 40 matched stars against 440 catalog stars, which is the replacement for the chunked-vs-unchunked check the chunking approach needed.

Tests/TestMatchStars.py, TestMath.py and TestSatellitePositions.py still pass alongside it.

1 and 4

Description updated for 4 — the units half is a false-positive reduction, not a no-op, and the 23/17 diff is mostly re-indentation.

Finding 1 (the magnitude pre-filter selecting mag < catalog_limit + 0.3, i.e. everything) is real and I've noted it in the description, but I'd rather not fold it in here: both call sites pass the catalog's own depth as lim_mag, so fixing it changes which stars get rejected as blended, and that deserves its own before/after on real station data rather than riding along with a units-and-memory fix. Will open it separately.

Not verified

No live station run — this is desk-verified plus the test suite. The blend counts should shift slightly with the wider cone on wide-lens stations (more catalog stars considered, so if anything marginally more blends caught), and I have no measurement of that on real data.

@Cybis320

Copy link
Copy Markdown
Contributor Author

Follow-up on the SkyFit2 half of the duplication point — 608a5c89 deletes that cone instead of sharing it, which is a better answer than the helper swap I pushed earlier.

filterCatalogStarsInsideFOV already takes a lim_mag kwarg on prerelease (SkyFit2.py:12299-12311), so the comment justifying the hand-rolled copy — "not self.filterCatalogStarsInsideFOV which incorrectly uses self.cat_lim_mag instead of the test catalog's LM" — is out of date. Passing lim_mag=test_lm gets the intended behaviour from the existing method:

_, test_catalog = self.filterCatalogStarsInsideFOV(test_catalog, lim_mag=test_lm)

That gives three things at once: it centres the cone via computeCentreRADec() at the current image time (your finding 2), it uses getFOVSelectionRadius internally, and it applies the test catalog's own LM rather than self.cat_lim_mag. Net −18 lines, and the shared helper is no longer needed at this call site — catalogStarsInFOV now has exactly one caller, filterBlendedStars.

The epochs match: jd in the tuning routine is bound at SkyFit2.py:5040 as date2JD(*self.img_handle.currentTime()), which is the same source computeCentreRADec() reads, so this is not a change of reference time — only of which pointing the cone is built around.

Worth flagging where this comes from, since it's stronger evidence than the synthetic test I added for finding 2. We hit this in production on our test-fleet branch: after a whole-night platepar refit, platepar.JD sat hours away from the displayed frame, the cone pointed at the wrong sky, the genuine FOV stars were all excluded, and the only survivors were back-projection folds — Phase 3 of the tuner matched essentially nothing against a perfectly good platepar. Fixed there by this same call. So finding 2 is not theoretical on an alt-az station; it just needed a large enough epoch gap to become visible.

Two related notes from merging this branch into that fleet branch, neither affecting this PR as it stands:

  • Our copy of filterBlendedStars has since grown a physical blend criterion — a neighbour only counts if its flux ratio at that separation pulls the photocentre more than ~0.5 px, since at coarse plate scales 2x FWHM is tens of arcminutes and some faint neighbour is nearly always inside it. If that lands upstream later, the cKDTree k-nearest query here has to become a radius query: the criterion depends on brightness, so a more distant but brighter neighbour can be the threatening one and k-nearest can miss it. Not an issue for the distance-only criterion on this branch.
  • That same criterion turns your finding 1 into something actionable — it yields a real magnitude cutoff (pull_px/(r_max - pull_px) as a flux-ratio floor, converted to a magnitude), which is a principled replacement for the lim_mag + mag_margin test that currently selects everything. That's the shape the separate PR should probably take.

@dvida

dvida commented Aug 6, 2026

Copy link
Copy Markdown
Contributor

Went through the final state carefully. One correctness regression to fix, plus a stale description; after those I'm happy to merge. Your follow-up comments hold up, including where they contradict the earlier review.

Blocking

B1 — k=3 nearest-neighbour can silently miss a blend

The old code scanned all catalog stars for 0.1 < d < blend_radius. The new code inspects only the 3 nearest and applies the same test. The d > 0.1 test is there to drop the star's own catalog entry, which projects to d = 0 exactly.

Failure mode: if three or more catalog rows sit within 0.1 px of a matched star, they fill all three k slots and a genuine neighbour at, say, 3 px is never examined — the blended star survives the filter. At the 4 mm plate scale 0.1 px is about 6.4 arcsec, and a deep Gaia catalog has plenty of triples inside 6.4" in dense fields. That's precisely the heavily-blended case the filter exists to catch, and it fails silently.

test_duplicate_catalog_entry_does_not_hide_a_blend covers two coincident rows (self + one duplicate), which is exactly the case k=3 survives. Three is not covered.

I'd swap the k-nearest query for a radius query. It is exactly equivalent to the old brute-force result, still bounded in memory, and it drops the k heuristic entirely. edge_margin is already this same radius, so the bounds cull and the search read as one idea rather than two:

        pts = np.column_stack([all_matched_x, all_matched_y])
        neighbour_lists = tree.query_ball_point(pts, edge_margin)

        # Neighbours within each star's own blend radius, excluding its own catalog entry
        #   (and any duplicate of it), which project to d ~ 0
        has_neighbor = np.zeros(len(check_indices), dtype=bool)
        for k, idx_list in enumerate(neighbour_lists):
            if not idx_list:
                continue
            d = np.hypot(catalog_x[idx_list] - all_matched_x[k],
                         catalog_y[idx_list] - all_matched_y[k])
            has_neighbor[k] = np.any((d < blend_radii[k]) & (d > 0.1))

Note: use the scalar r form. query_ball_point with an array r is a much later scipy addition and requirements.txt only pins scipy>=1.0.0.

This also removes the np.reshape(..., (len(check_indices), n_neighbors)) step, which only exists to normalise scipy's output shape when k == 1.

Worth adding a test alongside it: the matched star's position in the catalog three times plus a real neighbour at 4.5 px. It fails on the current code and passes after the change. test_matches_bruteforce_reference should pass unchanged — that's the one that proves the radius query is equivalent to the old matrix.

Incidentally this also removes the obstacle you flagged for the physical blend criterion on your fleet branch: a radius query is already the shape that criterion needs, since a more distant but brighter neighbour can be the threatening one and k-nearest can miss it.

B2 — Description is stale, and it becomes the squash message

The body still says the helper is "called from filterBlendedStars and from SkyFit2.count_matches_at_lm". 608a5c89 deleted that second caller in favour of filterCatalogStarsInsideFOV, so catalogStarsInFOV now has exactly one. The body also never describes the SkyFit2 rewrite at all, which is the more interesting half of the change. Please bring it in line with the final state.

Non-blocking

N1 — test_deep_catalog_outside_the_fov_is_cheap asserts something it doesn't construct

It scatters 20000 uniform whole-sky stars at mag 5.0 (all pass the magnitude cut) and asserts removed == 1. But the 4 mm template images roughly 3772 deg², about 9% of the sky, so ~1800 of those 20000 land inside the frame, at ~0.002 stars/px². With a 6 px blend radius that's ~0.22 expected neighbours per matched star, ~1.1 across the five — so there's roughly a 2-in-3 chance a different seed gives removed > 1. It passes on RandomState(7) by luck rather than by construction, and those in-frame stars are legitimate blends, not the "cannot possibly blend" stars the docstring describes.

Rejecting sampled stars within getFOVSelectionRadius(pp)*1.5 of the FOV centre before adding them makes removed == 1 guaranteed and makes the test mean what its name says.

N2 — Undocumented third behaviour change: the horizon cut

filterCatalogStarsInsideFOVsubsetCatalog(..., remove_under_horizon=True) drops catalog stars below −5 deg elevation. The hand-rolled cone it replaces never did. For a low-pointing camera whose FOV clips the horizon, n_catalog shrinks and the LM optimum in _findOptimalCatalogLM can shift. I think this is desirable — it makes the tuner agree with what SkyFit2 actually displays — but it's a third behaviour change beyond the two the description lists, so it should be named there.

N3 — min(fov_radius*margin, 90) can shrink the cone below the true FOV

On a fisheye platepar where getFOVSelectionRadius() exceeds 90 deg, the cap clips genuinely in-frame stars. Pre-existing — the old code capped at 90 too — so not a regression, but the new helper is the natural place to stop propagating it. Your call whether to fold it in here or leave it.

N4 — Cosmetic

Docstrings elsewhere in StarFilters.py put the summary on the line after """; the new helper puts it on the same line. Repo-wide the same-line form dominates, so either is fine, but pick one within the file.

Keep catalogStarsInFOV even though it's down to one caller — it's the unit the tests exercise and it names a real concept.

What I checked and agree with

  • The units bug and the fix are both real. getFOVSelectionRadius() projects the four image corners through the platepar with distortion and takes the max separation from centre, so it is exact by construction rather than an estimate.
  • Your correction on the arithmetic is the one that's right: the swap is not a size reduction. On the shipped 4 mm template the exact corner radius (55.7 deg) is wider than the linear estimate (47.0 deg), so the linear form would clip stars genuinely in frame. test_image_corners_are_inside at margin=1.001 is the right guard for that.
  • Centring the cone on jd is a bigger fix than the description claims. platepar.RA_d/dec_d is the raw reference pointing at platepar.JD; both cyXYToRADec and cyraDecToXY route through pointingCorrection() (CyFunctions.pyx:1298), which corrects that reference to the target JD. So the old code compared an uncorrected reference pointing against catalog coordinates, while the new xyToRaDecPP() centre sits in the same frame the catalog is projected from. Frame-consistent, not just time-consistent — which also explains why the production symptom you saw was as severe as it was.
  • Both filterBlendedStars call sites pass the image JD, so the jd-centred cone is correct at both: SkyFit2.py:8702 via currentFrameTime(), and AutoPlatepar.py:657.
  • The OOM path is the one being fixed. SkyFit2 passes catalog_stars_filtered_unmasked, already FOV-subset, so the cone is a no-op there — but AutoPlatepar.py:657 passes the full catalog. That's where the matrices blew up.
  • The filterCatalogStarsInsideFOV swap is sound, and I checked the preconditions that weren't discussed: subsetCatalog (CyFunctions.pyx:114) breaks on descending declination so it requires dec-sorted input, and readStarCatalog sorts descending (StarCatalog.py:519); its second mag <= test_lm cut is exactly redundant with readStarCatalog(lim_mag=test_lm) on the same synthetic-band column, so no stars are silently lost; and margin 1.0 there is correct rather than sloppy, since the next step is an in_image bounds check and the radius circumscribes the corners.
  • Deferring the magnitude pre-filter is the right call, and it's now safe to defer: memory is bounded by the cone, the bounds cull and the tree regardless of how many stars pass the magnitude test. The flux-ratio floor you describe sounds like the right shape for that PR.

Style matches the repo throughout — naming, docstring blocks, continuations, .format() in the library file, and the test file following Tests/TestMatchStars.py.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants