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