Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Cybis320 wants to merge 5 commits into
base: prerelease
Choose a base branch
from
Open

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

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
cd62e52
Fix the blend filter's FOV pre-filter units and bound its memory
Cybis320 Jul 24, 2026
1910164
Use the exact FOV selection radius and centre the catalog cone on the…
Cybis320 Jul 28, 2026
e41c4ad
Replace the blend distance matrix with an image-bounds cull and a KD-…
Cybis320 Jul 28, 2026
b812ff9
Add tests for the blended-star catalog pre-filters
Cybis320 Jul 28, 2026
608a5c8
Reuse filterCatalogStarsInsideFOV for the LM search catalog cone
Cybis320 Jul 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 20 additions & 11 deletions RMS/Astrometry/StarFilters.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

@dvida dvida Jul 28, 2026

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.

fov_radius = min(fov_radius, 90) # Cap at 90 degrees

in_fov = ang_dist_deg < fov_radius
Expand Down Expand Up @@ -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)

@dvida dvida Jul 28, 2026

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.

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)

@dvida dvida Jul 28, 2026

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.


# 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]:
Expand Down
3 changes: 2 additions & 1 deletion Utils/SkyFit2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

@dvida dvida Jul 28, 2026

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.

fov_radius = min(fov_radius, 90)

in_fov = ang_dist_deg < fov_radius
Expand Down