Fix the blend filter's FOV pre-filter units and bound its memory - #939
Fix the blend filter's FOV pre-filter units and bound its memory#939Cybis320 wants to merge 5 commits into
Conversation
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
left a comment
There was a problem hiding this comment.
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_magBoth call sites pass the catalog's own limiting magnitude as lim_mag:
RMS/Astrometry/AutoPlatepar.py:657-658passesconfig.catalog_mag_limit, and thecatalog_starsit forwards is the whole-sky catalog fromloadCatalogStars(config, config.catalog_mag_limit)(AutoPlatepar.py:784) — never FOV-subset.Utils/SkyFit2.py:8620-8628passesself.cat_lim_mag, which is the catalog load depth (SkyFit2.py:2483), and the LM search pushes it deep atSkyFit2.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_radius≥getFOVSelectionRadius(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+1against 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.
| # 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 |
There was a problem hiding this comment.
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.5Three reasons this is more than style:
- 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.
- It handles lenses the linear approximation only approximates.
F_scaleis the central scale; the corner angle depends on the distortion polynomial, whichgetFOVSelectionRadiusevaluates rather than assumes. - 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 callinggetFOVSelectionRadiusfrom both) makes the third copy impossible.
| # the process OOM-killed) | ||
| # Shape per chunk: (n_matched, chunk) | ||
| n_matched = len(check_indices) | ||
| chunk_size = max(1, int(5e6) // n_matched) |
There was a problem hiding this comment.
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
sqrtentirely: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.
| 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) |
There was a problem hiding this comment.
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.
| # 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 |
There was a problem hiding this comment.
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.
|
Thanks — took findings 2, 3 and 5 plus both inline suggestions, in The "halves the catalog again" number doesn't holdThe comparison is 76° with the 1.5 margin against ~50° without it, so it's measuring the margin, not the exactness. Keeping the Measured on the two platepars shipped in
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. 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
|
|
Follow-up on the SkyFit2 half of the duplication point —
_, test_catalog = self.filterCatalogStarsInsideFOV(test_catalog, lim_mag=test_lm)That gives three things at once: it centres the cone via The epochs match: 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, Two related notes from merging this branch into that fleet branch, neither affecting this PR as it stands:
|
|
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. BlockingB1 —
|
The FOV radius pre-filter in
filterBlendedStarsmultiplied the pixel diagonal byF_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()inRMS/Astrometry/StarFilters.py, which usesgetFOVSelectionRadius()— 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 ownjdrather than atplatepar.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
cKDTreenearest-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 aslim_mag, sobright_maskselectsmag < 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).