Skip to main content

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}