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