target_match/matcher.rs
1// This Source Code Form is subject to the terms of the Mozilla Public
2// License, v. 2.0. If a copy of the MPL was not distributed with this
3// file, You can obtain one at https://mozilla.org/MPL/2.0/.
4
5//! The matching engine: input trait, constraints, ranking, and a prebuilt index.
6//!
7//! A consumer's catalogue type implements [`SkyObject`] (position only — never a
8//! name). A [`Constraint`] combines a [`Membership`] shape (circular, rectangle,
9//! or rotated rectangle) with a [`Query`] mode. [`rank`] scans a slice; [`Matcher`]
10//! builds a declination-sorted index once and answers repeated queries, returning
11//! results identical to [`rank`].
12//!
13//! # Geometry
14//!
15//! Matching precesses the pointing to J2000 (via [`skymath::precess`]), then works
16//! in the local tangent frame about the pointing: an object's offset is decomposed
17//! from its great-circle separation and position angle (East of North, both from
18//! `skymath`) into East/North components, which are rotated into the camera frame
19//! for rectangular membership. The circumscribed circle pre-filters both rectangle
20//! tests, so the tangent decomposition is only evaluated for objects near the
21//! frame — never on the far side of the sky.
22
23use core::cmp::Ordering;
24
25use skymath::{position_angle, precess, separation, Angle, Epoch, Equatorial};
26
27use crate::optics::{Field, RadiusPolicy};
28
29/// A catalogue object that can be matched by sky position.
30///
31/// The trait exposes **only** a J2000 position — matching never reads a name or
32/// designation. A caller's own type keeps its identity; a [`Match`] borrows it.
33///
34/// # Example
35///
36/// ```
37/// use skymath::{Angle, Equatorial};
38/// use target_match::SkyObject;
39///
40/// struct Target {
41/// name: &'static str,
42/// ra: f64,
43/// dec: f64,
44/// }
45/// impl SkyObject for Target {
46/// fn position(&self) -> Equatorial {
47/// Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
48/// }
49/// }
50///
51/// let m31 = Target { name: "M 31", ra: 10.6847, dec: 41.2688 };
52/// assert!((m31.position().ra().degrees() - 10.6847).abs() < 1e-9);
53/// ```
54pub trait SkyObject {
55 /// The object's J2000 equatorial position — the only thing [`rank`],
56 /// [`is_framed`], and [`Matcher`] read from an implementer.
57 ///
58 /// # Example
59 ///
60 /// ```
61 /// use skymath::{Angle, Equatorial};
62 /// use target_match::SkyObject;
63 ///
64 /// struct Target { ra_deg: f64, dec_deg: f64 }
65 /// impl SkyObject for Target {
66 /// fn position(&self) -> Equatorial {
67 /// Equatorial::j2000(Angle::from_degrees(self.ra_deg), Angle::from_degrees(self.dec_deg)).unwrap()
68 /// }
69 /// }
70 ///
71 /// let m31 = Target { ra_deg: 10.6847, dec_deg: 41.2688 };
72 /// assert!((m31.position().ra().degrees() - 10.6847).abs() < 1e-9);
73 /// ```
74 fn position(&self) -> Equatorial;
75}
76
77/// The shape that decides whether an object is "in frame".
78///
79/// Usually built for you by a [`Constraint`] constructor; pass one directly to
80/// [`is_framed`] to test a single object.
81///
82/// # Example
83///
84/// ```
85/// use skymath::{Angle, Equatorial, ParseMode};
86/// use target_match::{is_framed, Membership, SkyObject};
87///
88/// struct Target {
89/// ra: f64,
90/// dec: f64,
91/// }
92/// impl SkyObject for Target {
93/// fn position(&self) -> Equatorial {
94/// Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
95/// }
96/// }
97///
98/// let m31 = Target { ra: 10.6847, dec: 41.2688 };
99/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
100///
101/// let shape = Membership::Circular { radius: Angle::from_degrees(1.0) };
102/// assert!(is_framed(pointing, &m31, shape).in_frame);
103/// ```
104#[derive(Debug, Clone, Copy, PartialEq)]
105#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
106pub enum Membership {
107 /// Pure angular distance: in frame iff separation ≤ `radius`.
108 Circular {
109 /// Search radius.
110 radius: Angle,
111 },
112 /// Axis-aligned rectangle of the given width×height field of view.
113 Rectangle {
114 /// `(width, height)` field of view.
115 fov: (Angle, Angle),
116 },
117 /// Rectangle rotated by a camera position angle (degrees East of North).
118 Rotated {
119 /// `(width, height)` field of view.
120 fov: (Angle, Angle),
121 /// Camera position angle, East of North.
122 position_angle: Angle,
123 },
124}
125
126/// What to return from a match.
127///
128/// Set on a [`Constraint`] via [`Constraint::all`], [`Constraint::nearest_one`],
129/// [`Constraint::nearest_n`], or [`Constraint::nearest_n_within`].
130///
131/// # Example
132///
133/// ```
134/// use skymath::{Angle, Equatorial, ParseMode};
135/// use target_match::{rank, Constraint, SkyObject};
136///
137/// struct Target {
138/// ra: f64,
139/// dec: f64,
140/// }
141/// impl SkyObject for Target {
142/// fn position(&self) -> Equatorial {
143/// Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
144/// }
145/// }
146///
147/// let catalog = [Target { ra: 10.6847, dec: 41.2688 }];
148/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
149///
150/// // `nearest_n(1)` sets the query to `Query::NearestN { n: 1, .. }`.
151/// let hits = rank(pointing, &catalog, Constraint::circular(Angle::from_degrees(1.0)).nearest_n(1));
152/// assert_eq!(hits.len(), 1);
153/// ```
154#[derive(Debug, Clone, Copy, PartialEq)]
155#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
156pub enum Query {
157 /// Every object inside the membership shape, ranked nearest-first.
158 AllWithinField,
159 /// The single nearest in-frame object.
160 NearestOne,
161 /// The `n` nearest objects by separation, optionally bounded by a radius.
162 NearestN {
163 /// Maximum number of results.
164 n: usize,
165 /// Optional maximum separation; unbounded when `None`.
166 max_radius: Option<Angle>,
167 },
168}
169
170/// A [`Membership`] shape combined with a [`Query`] mode (and the plate scale,
171/// when known, so pixel offsets can be reported).
172///
173/// Build one with a shape constructor ([`within`](Constraint::within),
174/// [`circular`](Constraint::circular), [`frame`](Constraint::frame),
175/// [`frame_rotated`](Constraint::frame_rotated)), then optionally switch the
176/// query mode ([`all`](Constraint::all), [`nearest_one`](Constraint::nearest_one),
177/// [`nearest_n`](Constraint::nearest_n),
178/// [`nearest_n_within`](Constraint::nearest_n_within)). Pass the result to
179/// [`rank`], [`Matcher::query`], or [`is_framed`].
180///
181/// # Example
182///
183/// ```
184/// use skymath::{Angle, Equatorial, ParseMode};
185/// use target_match::{rank, Constraint, Field, Optics, RadiusPolicy, SkyObject};
186///
187/// struct Target {
188/// ra: f64,
189/// dec: f64,
190/// }
191/// impl SkyObject for Target {
192/// fn position(&self) -> Equatorial {
193/// Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
194/// }
195/// }
196///
197/// let catalog = [Target { ra: 10.6847, dec: 41.2688 }];
198/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
199/// let field = Field::from_optics(Optics {
200/// focal_mm: 800.0,
201/// pixel_um: (3.76, 3.76),
202/// binning: (1, 1),
203/// pixels: (6248, 4176),
204/// })
205/// .unwrap();
206///
207/// let c = Constraint::within(&field, RadiusPolicy::Circumscribed).nearest_one();
208/// assert_eq!(rank(pointing, &catalog, c).len(), 1);
209/// ```
210#[derive(Debug, Clone, Copy, PartialEq)]
211#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
212pub struct Constraint {
213 /// The in-frame shape.
214 pub membership: Membership,
215 /// The query mode.
216 pub query: Query,
217 /// Per-axis plate scale (arcsec/px), used only to fill pixel offsets.
218 pub pixel_scale: Option<(f64, f64)>,
219}
220
221impl Constraint {
222 /// Circular membership with a radius derived from a field under `policy`.
223 ///
224 /// # Example
225 ///
226 /// ```
227 /// use target_match::{Constraint, Field, Optics, RadiusPolicy};
228 ///
229 /// let field = Field::from_optics(Optics {
230 /// focal_mm: 800.0,
231 /// pixel_um: (3.76, 3.76),
232 /// binning: (1, 1),
233 /// pixels: (6248, 4176),
234 /// })
235 /// .unwrap();
236 /// let c = Constraint::within(&field, RadiusPolicy::Circumscribed);
237 /// assert!(c.pixel_scale.is_some(), "the field's plate scale carries through");
238 /// ```
239 #[must_use]
240 pub fn within(field: &Field, policy: RadiusPolicy) -> Self {
241 Self {
242 membership: Membership::Circular {
243 radius: field.radius(policy),
244 },
245 query: Query::AllWithinField,
246 pixel_scale: field.pixel_scale(),
247 }
248 }
249 /// Circular membership with an explicit radius (no plate scale).
250 ///
251 /// # Example
252 ///
253 /// ```
254 /// use skymath::Angle;
255 /// use target_match::Constraint;
256 ///
257 /// let c = Constraint::circular(Angle::from_degrees(2.0));
258 /// assert!(c.pixel_scale.is_none(), "no `Field`, so no plate scale");
259 /// ```
260 #[must_use]
261 pub fn circular(radius: Angle) -> Self {
262 Self {
263 membership: Membership::Circular { radius },
264 query: Query::AllWithinField,
265 pixel_scale: None,
266 }
267 }
268 /// Axis-aligned rectangular membership from a field's width×height.
269 ///
270 /// # Example
271 ///
272 /// ```
273 /// use target_match::{Constraint, Field, Membership, Optics};
274 ///
275 /// let field = Field::from_optics(Optics {
276 /// focal_mm: 800.0,
277 /// pixel_um: (3.76, 3.76),
278 /// binning: (1, 1),
279 /// pixels: (6248, 4176),
280 /// })
281 /// .unwrap();
282 /// let c = Constraint::frame(&field);
283 /// assert!(matches!(c.membership, Membership::Rectangle { .. }));
284 /// ```
285 #[must_use]
286 pub fn frame(field: &Field) -> Self {
287 Self {
288 membership: Membership::Rectangle {
289 fov: (field.width(), field.height()),
290 },
291 query: Query::AllWithinField,
292 pixel_scale: field.pixel_scale(),
293 }
294 }
295 /// Rotated rectangular membership from a field and a camera position angle
296 /// (degrees East of North).
297 ///
298 /// # Example
299 ///
300 /// An object due north of the pointing sits on the frame's `+y` axis when
301 /// axis-aligned; a 90° East-of-North camera rotation moves it onto the `-x`
302 /// axis.
303 ///
304 /// ```
305 /// use skymath::{Angle, Equatorial, ParseMode};
306 /// use target_match::{rank, Constraint, Field, SkyObject};
307 ///
308 /// struct Target {
309 /// ra: f64,
310 /// dec: f64,
311 /// }
312 /// impl SkyObject for Target {
313 /// fn position(&self) -> Equatorial {
314 /// Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
315 /// }
316 /// }
317 ///
318 /// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
319 /// let field = Field::from_fov(Angle::from_degrees(2.0), Angle::from_degrees(2.0)).unwrap();
320 /// let catalog = [Target { ra: 10.6847, dec: 41.2688 + 0.3 }]; // ~0.3° due north
321 ///
322 /// let (ax, ay) = rank(pointing, &catalog, Constraint::frame(&field))[0].offset.frame.unwrap();
323 /// assert!(ay.degrees() > 0.25 && ax.degrees().abs() < 0.05);
324 ///
325 /// let rotated = Constraint::frame_rotated(&field, Angle::from_degrees(90.0));
326 /// let (rx, ry) = rank(pointing, &catalog, rotated)[0].offset.frame.unwrap();
327 /// assert!(rx.degrees() < -0.25 && ry.degrees().abs() < 0.05);
328 /// ```
329 #[must_use]
330 pub fn frame_rotated(field: &Field, position_angle: Angle) -> Self {
331 Self {
332 membership: Membership::Rotated {
333 fov: (field.width(), field.height()),
334 position_angle,
335 },
336 query: Query::AllWithinField,
337 pixel_scale: field.pixel_scale(),
338 }
339 }
340 /// Set the query to all-within-field.
341 ///
342 /// # Example
343 ///
344 /// ```
345 /// use skymath::Angle;
346 /// use target_match::{Constraint, Query};
347 ///
348 /// let c = Constraint::circular(Angle::from_degrees(1.0)).nearest_one().all();
349 /// assert_eq!(c.query, Query::AllWithinField);
350 /// ```
351 #[must_use]
352 pub fn all(mut self) -> Self {
353 self.query = Query::AllWithinField;
354 self
355 }
356 /// Set the query to nearest-one.
357 ///
358 /// # Example
359 ///
360 /// ```
361 /// use skymath::Angle;
362 /// use target_match::{Constraint, Query};
363 ///
364 /// let c = Constraint::circular(Angle::from_degrees(1.0)).nearest_one();
365 /// assert_eq!(c.query, Query::NearestOne);
366 /// ```
367 #[must_use]
368 pub fn nearest_one(mut self) -> Self {
369 self.query = Query::NearestOne;
370 self
371 }
372 /// Set the query to the `n` nearest (unbounded).
373 ///
374 /// # Example
375 ///
376 /// ```
377 /// use skymath::Angle;
378 /// use target_match::{Constraint, Query};
379 ///
380 /// let c = Constraint::circular(Angle::from_degrees(1.0)).nearest_n(3);
381 /// assert_eq!(c.query, Query::NearestN { n: 3, max_radius: None });
382 /// ```
383 #[must_use]
384 pub fn nearest_n(mut self, n: usize) -> Self {
385 self.query = Query::NearestN {
386 n,
387 max_radius: None,
388 };
389 self
390 }
391 /// Set the query to the `n` nearest within `max_radius`.
392 ///
393 /// # Example
394 ///
395 /// ```
396 /// use skymath::Angle;
397 /// use target_match::{Constraint, Query};
398 ///
399 /// let radius = Angle::from_degrees(1.0);
400 /// let c = Constraint::circular(radius).nearest_n_within(3, radius);
401 /// assert_eq!(c.query, Query::NearestN { n: 3, max_radius: Some(radius) });
402 /// ```
403 #[must_use]
404 pub fn nearest_n_within(mut self, n: usize, max_radius: Angle) -> Self {
405 self.query = Query::NearestN {
406 n,
407 max_radius: Some(max_radius),
408 };
409 self
410 }
411}
412
413/// The offset of a matched object relative to the frame centre.
414///
415/// Carried on every [`Match`]. `frame` and `pixels` are only populated for
416/// rectangular [`Membership`] (see [`Constraint::frame`],
417/// [`Constraint::frame_rotated`]) — circular membership has no frame axes.
418///
419/// # Example
420///
421/// ```
422/// use skymath::{Angle, Equatorial, ParseMode};
423/// use target_match::{rank, Constraint, Field, Optics, SkyObject};
424///
425/// struct Target {
426/// ra: f64,
427/// dec: f64,
428/// }
429/// impl SkyObject for Target {
430/// fn position(&self) -> Equatorial {
431/// Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
432/// }
433/// }
434///
435/// let catalog = [Target { ra: 10.6847, dec: 41.2688 }];
436/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
437/// let field = Field::from_optics(Optics {
438/// focal_mm: 800.0,
439/// pixel_um: (3.76, 3.76),
440/// binning: (1, 1),
441/// pixels: (6248, 4176),
442/// })
443/// .unwrap();
444///
445/// let hits = rank(pointing, &catalog, Constraint::frame(&field));
446/// let offset = hits[0].offset;
447/// assert!(offset.frame.is_some(), "rectangular membership reports a frame-aligned offset");
448/// assert!(offset.pixels.is_some(), "plate scale from `Optics` fills the pixel offset");
449/// ```
450#[derive(Debug, Clone, Copy, PartialEq)]
451#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
452pub struct Offset {
453 /// Sky-tangent offset `(East, North)` — always present.
454 pub sky: (Angle, Angle),
455 /// Frame-aligned offset `(x, y)`, present for rectangular membership.
456 pub frame: Option<(Angle, Angle)>,
457 /// Frame-aligned offset in pixels `(x, y)`, present when a plate scale is known.
458 pub pixels: Option<(f64, f64)>,
459}
460
461/// A ranked match: a borrowed catalogue object plus its computed geometry.
462///
463/// Returned by [`rank`], [`Matcher::query`], and [`is_framed`].
464///
465/// # Example
466///
467/// ```
468/// use skymath::{Angle, Equatorial, ParseMode};
469/// use target_match::{rank, Constraint, Field, Optics, RadiusPolicy, SkyObject};
470///
471/// struct Target {
472/// name: &'static str,
473/// ra: f64,
474/// dec: f64,
475/// }
476/// impl SkyObject for Target {
477/// fn position(&self) -> Equatorial {
478/// Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
479/// }
480/// }
481///
482/// let catalog = [Target { name: "M 31", ra: 10.6847, dec: 41.2688 }];
483/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
484/// let field = Field::from_optics(Optics {
485/// focal_mm: 800.0,
486/// pixel_um: (3.76, 3.76),
487/// binning: (1, 1),
488/// pixels: (6248, 4176),
489/// })
490/// .unwrap();
491///
492/// let hits = rank(pointing, &catalog, Constraint::within(&field, RadiusPolicy::Circumscribed));
493/// let m = &hits[0];
494/// assert_eq!(m.object.name, "M 31");
495/// assert!(m.in_frame);
496/// assert!(m.separation.arcseconds() < 2.0);
497/// ```
498#[derive(Debug, Clone, Copy, PartialEq)]
499pub struct Match<'a, T> {
500 /// The matched object, borrowed from the caller's slice or the [`Matcher`].
501 pub object: &'a T,
502 /// Great-circle separation from the frame centre.
503 pub separation: Angle,
504 /// Whether the object is inside the active membership shape.
505 pub in_frame: bool,
506 /// The object's offset relative to the frame centre.
507 pub offset: Offset,
508 /// Position angle from frame centre to the object (degrees East of North).
509 pub position_angle: Angle,
510}
511
512/// Tangent-frame `(East, North)` offset in radians, decomposed from an
513/// already-computed separation and position angle (the polar decomposition
514/// `skymath::tangent_offset` performs, reusing this evaluation's sep/PA
515/// instead of recomputing them).
516fn tangent_components(sep: Angle, pa: Angle) -> (f64, f64) {
517 let (s, r) = (pa.radians(), sep.radians());
518 (r * s.sin(), r * s.cos())
519}
520
521fn rotate(east: f64, north: f64, pa_rad: f64) -> (f64, f64) {
522 let (s, c) = pa_rad.sin_cos();
523 (east * c - north * s, east * s + north * c)
524}
525
526/// The circumscribed-circle radius (radians) that bounds a membership shape.
527fn bound_radius(m: Membership) -> f64 {
528 match m {
529 Membership::Circular { radius } => radius.radians(),
530 Membership::Rectangle { fov } | Membership::Rotated { fov, .. } => {
531 (fov.0.radians() / 2.0).hypot(fov.1.radians() / 2.0)
532 }
533 }
534}
535
536fn contains(sep: Angle, east: f64, north: f64, m: Membership) -> bool {
537 match m {
538 Membership::Circular { radius } => {
539 let r = radius.radians();
540 r.is_finite() && r >= 0.0 && sep.radians() <= r
541 }
542 Membership::Rectangle { fov } => within_rect(sep, east, north, fov, 0.0),
543 Membership::Rotated {
544 fov,
545 position_angle,
546 } => within_rect(sep, east, north, fov, position_angle.radians()),
547 }
548}
549
550fn within_rect(sep: Angle, east: f64, north: f64, fov: (Angle, Angle), pa: f64) -> bool {
551 let (hx, hy) = (fov.0.radians() / 2.0, fov.1.radians() / 2.0);
552 let circum = hx.hypot(hy);
553 // Circumscribed pre-filter (also rejects NaN separations).
554 if sep.radians() > circum || sep.radians().is_nan() {
555 return false;
556 }
557 let (x, y) = rotate(east, north, pa);
558 x.abs() <= hx && y.abs() <= hy
559}
560
561fn build_offset(east: f64, north: f64, m: Membership, scale: Option<(f64, f64)>) -> Offset {
562 let sky = (Angle::from_radians(east), Angle::from_radians(north));
563 match m {
564 Membership::Circular { .. } => Offset {
565 sky,
566 frame: None,
567 pixels: None,
568 },
569 Membership::Rectangle { .. } | Membership::Rotated { .. } => {
570 let pa = match m {
571 Membership::Rotated { position_angle, .. } => position_angle.radians(),
572 _ => 0.0,
573 };
574 let (x, y) = rotate(east, north, pa);
575 let frame = Some((Angle::from_radians(x), Angle::from_radians(y)));
576 let pixels = scale.map(|(sx, sy)| {
577 (
578 Angle::from_radians(x).arcseconds() / sx,
579 Angle::from_radians(y).arcseconds() / sy,
580 )
581 });
582 Offset { sky, frame, pixels }
583 }
584 }
585}
586
587fn evaluate<'a, T: SkyObject>(
588 pointing: Equatorial,
589 obj: &'a T,
590 m: Membership,
591 scale: Option<(f64, f64)>,
592) -> Match<'a, T> {
593 let pos = obj.position();
594 let sep = separation(pointing, pos);
595 let pa = position_angle(pointing, pos);
596 let (east, north) = tangent_components(sep, pa);
597 Match {
598 object: obj,
599 separation: sep,
600 in_frame: contains(sep, east, north, m),
601 offset: build_offset(east, north, m, scale),
602 position_angle: pa,
603 }
604}
605
606/// Whether `obj` should be kept for `query` given its evaluated match.
607fn keep<T>(m: &Match<'_, T>, query: Query) -> bool {
608 match query {
609 Query::AllWithinField | Query::NearestOne => m.in_frame,
610 Query::NearestN { max_radius, .. } => {
611 max_radius.map_or(true, |r| m.separation.radians() <= r.radians())
612 }
613 }
614}
615
616/// Shared core: evaluate candidates, filter, rank, and truncate per the query.
617fn rank_candidates<'a, T: SkyObject, I>(
618 pointing: Equatorial,
619 candidates: I,
620 c: &Constraint,
621) -> Vec<Match<'a, T>>
622where
623 I: Iterator<Item = (usize, &'a T)>,
624{
625 let mut scored: Vec<(usize, Match<'a, T>)> = candidates
626 .map(|(i, o)| (i, evaluate(pointing, o, c.membership, c.pixel_scale)))
627 .filter(|(_, m)| keep(m, c.query))
628 .collect();
629 scored.sort_by(|a, b| {
630 a.1.separation
631 .radians()
632 .partial_cmp(&b.1.separation.radians())
633 .unwrap_or(Ordering::Equal)
634 .then(a.0.cmp(&b.0))
635 });
636 let mut out: Vec<Match<'a, T>> = scored.into_iter().map(|(_, m)| m).collect();
637 match c.query {
638 Query::NearestOne => out.truncate(1),
639 Query::NearestN { n, .. } => out.truncate(n),
640 Query::AllWithinField => {}
641 }
642 out
643}
644
645/// Rank a slice of objects against a pointing under a constraint (stateless scan).
646///
647/// The pointing is precessed to J2000 first. Results are ascending by separation
648/// with ties broken by input order. For repeated queries against one catalogue,
649/// build a [`Matcher`] instead — it returns identical results faster.
650///
651/// # Example
652///
653/// ```
654/// use skymath::{Angle, Equatorial, ParseMode};
655/// use target_match::{rank, Constraint, Field, Optics, RadiusPolicy, SkyObject};
656///
657/// struct Target {
658/// name: &'static str,
659/// ra: f64,
660/// dec: f64,
661/// }
662/// impl SkyObject for Target {
663/// fn position(&self) -> Equatorial {
664/// Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
665/// }
666/// }
667///
668/// let catalog = [
669/// Target { name: "M 31", ra: 10.6847, dec: 41.2688 },
670/// Target { name: "M 33", ra: 23.4621, dec: 30.6599 },
671/// ];
672/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
673/// let field = Field::from_optics(Optics {
674/// focal_mm: 800.0,
675/// pixel_um: (3.76, 3.76),
676/// binning: (1, 1),
677/// pixels: (6248, 4176),
678/// })
679/// .unwrap();
680///
681/// let hits = rank(pointing, &catalog, Constraint::within(&field, RadiusPolicy::Circumscribed).nearest_one());
682/// assert_eq!(hits[0].object.name, "M 31");
683/// ```
684#[must_use]
685pub fn rank<T: SkyObject>(pointing: Equatorial, objects: &[T], c: Constraint) -> Vec<Match<'_, T>> {
686 let p = precess(pointing, Epoch::J2000);
687 rank_candidates(p, objects.iter().enumerate(), &c)
688}
689
690/// Evaluate a single object against a frame, returning its membership + geometry.
691///
692/// The pointing is precessed to J2000 first. Unlike [`rank`], no filtering or
693/// ranking is applied — the returned [`Match`] always describes `object`.
694///
695/// # Example
696///
697/// ```
698/// use skymath::{Angle, Equatorial, ParseMode};
699/// use target_match::{is_framed, Membership, SkyObject};
700///
701/// struct Target {
702/// ra: f64,
703/// dec: f64,
704/// }
705/// impl SkyObject for Target {
706/// fn position(&self) -> Equatorial {
707/// Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
708/// }
709/// }
710///
711/// let m31 = Target { ra: 10.6847, dec: 41.2688 };
712/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
713///
714/// let m = is_framed(pointing, &m31, Membership::Circular { radius: Angle::from_degrees(1.0) });
715/// assert!(m.in_frame);
716/// ```
717#[must_use]
718pub fn is_framed<T: SkyObject>(
719 pointing: Equatorial,
720 object: &T,
721 membership: Membership,
722) -> Match<'_, T> {
723 let p = precess(pointing, Epoch::J2000);
724 evaluate(p, object, membership, None)
725}
726
727/// A prebuilt, declination-sorted index for repeated queries against one catalogue.
728///
729/// Produces results identical to [`rank`] for the same objects, pointing, and
730/// constraint — the index is a performance optimization only. Build it once, then
731/// [`query`](Matcher::query) many pointings.
732///
733/// # Example
734///
735/// ```
736/// use skymath::{Angle, Equatorial, ParseMode};
737/// use target_match::{Constraint, Field, Matcher, Membership, Optics, RadiusPolicy, SkyObject};
738///
739/// struct Target {
740/// name: &'static str,
741/// ra: f64,
742/// dec: f64,
743/// }
744/// impl SkyObject for Target {
745/// fn position(&self) -> Equatorial {
746/// Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
747/// }
748/// }
749///
750/// let matcher = Matcher::from_objects(vec![
751/// Target { name: "M 31", ra: 10.6847, dec: 41.2688 },
752/// Target { name: "M 33", ra: 23.4621, dec: 30.6599 },
753/// ]);
754/// assert_eq!(matcher.objects().len(), 2);
755///
756/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
757/// let field = Field::from_optics(Optics {
758/// focal_mm: 800.0,
759/// pixel_um: (3.76, 3.76),
760/// binning: (1, 1),
761/// pixels: (6248, 4176),
762/// })
763/// .unwrap();
764///
765/// let hits = matcher.query(pointing, Constraint::within(&field, RadiusPolicy::Circumscribed).nearest_one());
766/// assert_eq!(hits[0].object.name, "M 31");
767///
768/// let m = matcher.is_framed(
769/// pointing,
770/// &matcher.objects()[0],
771/// Membership::Circular { radius: Angle::from_degrees(1.0) },
772/// );
773/// assert!(m.in_frame);
774/// ```
775pub struct Matcher<T> {
776 storage: Vec<T>,
777 /// `(declination_degrees, original_index)` sorted by declination.
778 sorted: Vec<(f64, usize)>,
779}
780
781impl<T: SkyObject> Matcher<T> {
782 /// Build an index from a set of objects (original order is preserved for
783 /// tie-breaking and [`objects`](Matcher::objects)).
784 ///
785 /// # Example
786 ///
787 /// ```
788 /// use skymath::{Angle, Equatorial};
789 /// use target_match::{Matcher, SkyObject};
790 ///
791 /// struct Target { ra_deg: f64, dec_deg: f64 }
792 /// impl SkyObject for Target {
793 /// fn position(&self) -> Equatorial {
794 /// Equatorial::j2000(Angle::from_degrees(self.ra_deg), Angle::from_degrees(self.dec_deg)).unwrap()
795 /// }
796 /// }
797 ///
798 /// let matcher = Matcher::from_objects(vec![
799 /// Target { ra_deg: 10.6847, dec_deg: 41.2688 },
800 /// Target { ra_deg: 23.4621, dec_deg: 30.6599 },
801 /// ]);
802 /// assert_eq!(matcher.objects().len(), 2);
803 /// ```
804 #[must_use]
805 pub fn from_objects(objects: Vec<T>) -> Self {
806 let mut sorted: Vec<(f64, usize)> = objects
807 .iter()
808 .enumerate()
809 .map(|(i, o)| (o.position().dec().degrees(), i))
810 .collect();
811 sorted.sort_by(|a, b| {
812 a.0.partial_cmp(&b.0)
813 .unwrap_or(Ordering::Equal)
814 .then(a.1.cmp(&b.1))
815 });
816 Self {
817 storage: objects,
818 sorted,
819 }
820 }
821
822 /// The stored objects, in their original insertion order.
823 ///
824 /// # Example
825 ///
826 /// ```
827 /// use skymath::{Angle, Equatorial};
828 /// use target_match::{Matcher, SkyObject};
829 ///
830 /// struct Target { name: &'static str, ra_deg: f64, dec_deg: f64 }
831 /// impl SkyObject for Target {
832 /// fn position(&self) -> Equatorial {
833 /// Equatorial::j2000(Angle::from_degrees(self.ra_deg), Angle::from_degrees(self.dec_deg)).unwrap()
834 /// }
835 /// }
836 ///
837 /// let matcher = Matcher::from_objects(vec![Target { name: "M 31", ra_deg: 10.6847, dec_deg: 41.2688 }]);
838 /// assert_eq!(matcher.objects()[0].name, "M 31");
839 /// ```
840 #[must_use]
841 pub fn objects(&self) -> &[T] {
842 &self.storage
843 }
844
845 /// Query the index for a pointing under a constraint. Results are
846 /// identical to calling [`rank`] with [`objects`](Matcher::objects).
847 ///
848 /// # Example
849 ///
850 /// ```
851 /// use skymath::{Angle, Equatorial, ParseMode};
852 /// use target_match::{Constraint, Field, Matcher, Optics, RadiusPolicy, SkyObject};
853 ///
854 /// struct Target { name: &'static str, ra_deg: f64, dec_deg: f64 }
855 /// impl SkyObject for Target {
856 /// fn position(&self) -> Equatorial {
857 /// Equatorial::j2000(Angle::from_degrees(self.ra_deg), Angle::from_degrees(self.dec_deg)).unwrap()
858 /// }
859 /// }
860 ///
861 /// let matcher = Matcher::from_objects(vec![Target { name: "M 31", ra_deg: 10.6847, dec_deg: 41.2688 }]);
862 /// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
863 /// let field = Field::from_optics(Optics {
864 /// focal_mm: 800.0, pixel_um: (3.76, 3.76), binning: (1, 1), pixels: (6248, 4176),
865 /// })
866 /// .unwrap();
867 ///
868 /// let hits = matcher.query(pointing, Constraint::within(&field, RadiusPolicy::Circumscribed).nearest_one());
869 /// assert_eq!(hits[0].object.name, "M 31");
870 /// ```
871 #[must_use]
872 pub fn query(&self, pointing: Equatorial, c: Constraint) -> Vec<Match<'_, T>> {
873 let p = precess(pointing, Epoch::J2000);
874 let r = match c.query {
875 Query::NearestN { max_radius, .. } => max_radius.map_or(f64::INFINITY, |a| a.radians()),
876 _ => bound_radius(c.membership),
877 };
878 let idxs = self.band(p.dec().degrees(), r);
879 rank_candidates(p, idxs.into_iter().map(|i| (i, &self.storage[i])), &c)
880 }
881
882 /// Evaluate a single stored-or-external object against a frame (see [`is_framed`]).
883 ///
884 /// # Example
885 ///
886 /// ```
887 /// use skymath::{Angle, Equatorial, ParseMode};
888 /// use target_match::{Matcher, Membership, SkyObject};
889 ///
890 /// struct Target { ra_deg: f64, dec_deg: f64 }
891 /// impl SkyObject for Target {
892 /// fn position(&self) -> Equatorial {
893 /// Equatorial::j2000(Angle::from_degrees(self.ra_deg), Angle::from_degrees(self.dec_deg)).unwrap()
894 /// }
895 /// }
896 ///
897 /// let matcher = Matcher::from_objects(vec![Target { ra_deg: 10.6847, dec_deg: 41.2688 }]);
898 /// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
899 ///
900 /// let obj = &matcher.objects()[0];
901 /// let m = matcher.is_framed(pointing, obj, Membership::Circular { radius: Angle::from_degrees(1.0) });
902 /// assert!(m.in_frame);
903 /// ```
904 #[must_use]
905 pub fn is_framed<'a>(
906 &self,
907 pointing: Equatorial,
908 object: &'a T,
909 m: Membership,
910 ) -> Match<'a, T> {
911 is_framed(pointing, object, m)
912 }
913
914 /// Original indices whose declination lies within `r_rad` of `dec0_deg`.
915 fn band(&self, dec0_deg: f64, r_rad: f64) -> Vec<usize> {
916 if r_rad.is_infinite() && r_rad > 0.0 {
917 return self.sorted.iter().map(|&(_, i)| i).collect();
918 }
919 if !r_rad.is_finite() || r_rad < 0.0 {
920 return Vec::new();
921 }
922 let r_deg = r_rad.to_degrees();
923 let (lo, hi) = (dec0_deg - r_deg, dec0_deg + r_deg);
924 let start = self.sorted.partition_point(|&(d, _)| d < lo);
925 let end = self.sorted.partition_point(|&(d, _)| d <= hi);
926 self.sorted[start..end].iter().map(|&(_, i)| i).collect()
927 }
928}
929
930#[cfg(test)]
931mod tests {
932 use super::*;
933
934 #[derive(Clone)]
935 struct Obj {
936 name: &'static str,
937 ra: f64,
938 dec: f64,
939 }
940 impl SkyObject for Obj {
941 fn position(&self) -> Equatorial {
942 Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
943 }
944 }
945
946 fn catalog() -> Vec<Obj> {
947 vec![
948 Obj {
949 name: "M 31",
950 ra: 10.6847,
951 dec: 41.2688,
952 },
953 Obj {
954 name: "M 110",
955 ra: 10.0921,
956 dec: 41.6853,
957 },
958 Obj {
959 name: "M 33",
960 ra: 23.4621,
961 dec: 30.6599,
962 },
963 Obj {
964 name: "M 42",
965 ra: 83.8221,
966 dec: -5.3911,
967 },
968 ]
969 }
970
971 fn m31() -> Equatorial {
972 Equatorial::j2000(Angle::from_degrees(10.6847), Angle::from_degrees(41.2688)).unwrap()
973 }
974
975 #[test]
976 fn nearest_one_circular() {
977 let cat = catalog();
978 let c = Constraint::circular(Angle::from_degrees(2.0)).nearest_one();
979 let hits = rank(m31(), &cat, c);
980 assert_eq!(hits.len(), 1);
981 assert_eq!(hits[0].object.name, "M 31");
982 assert!(hits[0].separation.arcseconds() < 1.0);
983 }
984
985 #[test]
986 fn all_within_field_ranked() {
987 let cat = catalog();
988 // ~1° radius keeps M31 (self) + M110 (~0.62°); excludes M33/M42.
989 let c = Constraint::circular(Angle::from_degrees(1.0)).all();
990 let hits = rank(m31(), &cat, c);
991 assert_eq!(
992 hits.len(),
993 2,
994 "{hits:?}",
995 hits = hits.iter().map(|h| h.object.name).collect::<Vec<_>>()
996 );
997 assert_eq!(hits[0].object.name, "M 31");
998 assert_eq!(hits[1].object.name, "M 110");
999 assert!(hits[0].separation.radians() <= hits[1].separation.radians());
1000 }
1001
1002 #[test]
1003 fn nearest_n_bounds_and_counts() {
1004 let cat = catalog();
1005 let c = Constraint::circular(Angle::from_degrees(1.0)).nearest_n(3);
1006 let hits = rank(m31(), &cat, c);
1007 assert_eq!(hits.len(), 3, "top-3 by separation regardless of frame");
1008 assert_eq!(hits[0].object.name, "M 31");
1009 // Bounded nearest-N respects the radius.
1010 let c2 = Constraint::circular(Angle::from_degrees(1.0))
1011 .nearest_n_within(3, Angle::from_degrees(1.0));
1012 assert_eq!(rank(m31(), &cat, c2).len(), 2, "only M31 + M110 within 1°");
1013 }
1014
1015 #[test]
1016 fn coordinates_only_never_name() {
1017 // A far object literally named "M 31" must NOT match near M31's pointing.
1018 let cat = vec![
1019 Obj {
1020 name: "M 31",
1021 ra: 200.0,
1022 dec: -40.0,
1023 },
1024 Obj {
1025 name: "Some Galaxy",
1026 ra: 10.6847,
1027 dec: 41.2688,
1028 },
1029 ];
1030 let c = Constraint::circular(Angle::from_degrees(2.0)).nearest_one();
1031 let hits = rank(m31(), &cat, c);
1032 assert_eq!(hits.len(), 1);
1033 assert_eq!(hits[0].object.name, "Some Galaxy");
1034 }
1035
1036 #[test]
1037 fn rectangle_excludes_circle_only_corner() {
1038 // An object just outside the rectangle but inside the circumscribed circle.
1039 // Field 2°×1°: half-width 1°, half-height 0.5°. Put an object 0.9° North
1040 // (in frame) vs 0.9° along the diagonal (out of the axis-aligned rect).
1041 let field = Field::from_fov(Angle::from_degrees(2.0), Angle::from_degrees(1.0)).unwrap();
1042 let north_obj = Obj {
1043 name: "N",
1044 ra: 10.6847,
1045 dec: 41.2688 + 0.4,
1046 }; // within 0.5° height
1047 let high_obj = Obj {
1048 name: "H",
1049 ra: 10.6847,
1050 dec: 41.2688 + 0.9,
1051 }; // beyond height, inside circle
1052 let cat = vec![north_obj, high_obj];
1053 let c = Constraint::frame(&field).all();
1054 let hits = rank(m31(), &cat, c);
1055 assert_eq!(hits.len(), 1);
1056 assert_eq!(hits[0].object.name, "N");
1057 }
1058
1059 #[test]
1060 fn matcher_matches_rank_exactly() {
1061 let cat = catalog();
1062 let c = Constraint::circular(Angle::from_degrees(5.0)).all();
1063 let via_rank: Vec<_> = rank(m31(), &cat, c).iter().map(|m| m.object.name).collect();
1064 let matcher = Matcher::from_objects(cat.clone());
1065 let via_index: Vec<_> = matcher
1066 .query(m31(), c)
1067 .iter()
1068 .map(|m| m.object.name)
1069 .collect();
1070 assert_eq!(via_rank, via_index);
1071 assert_eq!(matcher.objects().len(), 4);
1072 }
1073
1074 #[test]
1075 fn is_framed_reports_geometry() {
1076 let m110 = catalog()[1].clone();
1077 let m = is_framed(
1078 m31(),
1079 &m110,
1080 Membership::Circular {
1081 radius: Angle::from_degrees(1.0),
1082 },
1083 );
1084 assert!(m.in_frame);
1085 assert!((0.4..0.9).contains(&m.separation.degrees()));
1086 }
1087
1088 #[test]
1089 fn empty_catalog_and_zero_radius() {
1090 let cat = catalog();
1091 assert!(rank(
1092 m31(),
1093 &[] as &[Obj],
1094 Constraint::circular(Angle::from_degrees(1.0))
1095 )
1096 .is_empty());
1097 let c = Constraint::circular(Angle::from_degrees(-1.0)).all();
1098 assert!(
1099 rank(m31(), &cat, c).is_empty(),
1100 "negative radius matches nothing"
1101 );
1102 }
1103}