tetra3 0.1.0

Rust implementation of Tetra3: Fast and robust star plate solver
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230

//! Star plate solver based on the tetra3/cedar-solve algorithm.
//!
//! This module implements a "lost-in-space" plate solver that identifies star patterns
//! in images and determines the camera pointing direction. The algorithm:
//!
//! 1. **Database generation**: Precomputes 4-star geometric patterns from a star catalog,
//!    hashing their edge ratios into a lookup table.
//! 2. **Solving**: Given image centroids and an approximate FOV, tries 4-centroid
//!    combinations, looks up matching catalog patterns, estimates rotation via SVD,
//!    and verifies by counting star matches.
//!
//! Reference: cedar-solve / tetra3 by Gustav Pettersson (ESA) and Steven Rosenthal.

pub mod combinations;
pub mod database;
pub mod pattern;
pub mod solve;

use rkyv::{Archive, Deserialize, Serialize};

use crate::{Quaternion, StarCatalog};

// ── Status codes (matching tetra3) ──────────────────────────────────────────

/// Outcome of a plate-solve attempt.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum SolveStatus {
    /// A valid match was found.
    MatchFound,
    /// All pattern combinations were exhausted without a match.
    NoMatch,
    /// The solve timeout was reached before a match was found.
    Timeout,
    /// Too few centroids were provided to form a pattern.
    TooFew,
}

// ── Database properties ─────────────────────────────────────────────────────

/// Metadata describing how a solver database was built.
#[derive(Debug, Clone, Archive, Serialize, Deserialize)]
pub struct DatabaseProperties {
    /// Number of quantization bins per edge-ratio dimension.
    /// Computed as round(0.25 / pattern_max_error).
    pub pattern_bins: u32,
    /// Maximum tolerated error in edge ratios for a match.
    pub pattern_max_error: f32,
    /// Maximum FOV the database was built for (radians).
    pub max_fov_rad: f32,
    /// Minimum FOV the database was built for (radians).
    pub min_fov_rad: f32,
    /// Faintest star magnitude included.
    pub star_max_magnitude: f32,
    /// Total number of distinct patterns stored.
    pub num_patterns: u32,
    /// Catalog coordinate epoch (e.g. 2000 for J2000/ICRS).
    pub epoch_equinox: u16,
    /// Year to which proper motions were propagated.
    pub epoch_proper_motion_year: f32,
    /// Max catalog stars per FOV used for verification.
    pub verification_stars_per_fov: u32,
    /// Lattice field oversampling factor during generation.
    pub lattice_field_oversampling: u32,
    /// Number of patterns generated per lattice field.
    pub patterns_per_lattice_field: u32,
}

// ── The solver database ─────────────────────────────────────────────────────

/// Complete solver database, serializable with rkyv.
///
/// Contains a spatial star catalog (brightness-sorted), precomputed unit vectors,
/// and a pattern hash table for fast geometric matching.
#[derive(Debug, Clone, Archive, Serialize, Deserialize)]
pub struct SolverDatabase {
    /// Spatial star catalog. Stars are sorted brightest-first so that
    /// index order equals brightness order.
    pub star_catalog: StarCatalog,

    /// Precomputed ICRS unit vectors for each star, matching star_catalog.stars ordering.
    /// Stored as [x, y, z] where x=cos(ra)cos(dec), y=sin(ra)cos(dec), z=sin(dec).
    pub star_vectors: Vec<[f32; 3]>,

    /// Original catalog IDs (e.g. HIP number) for each star.
    pub star_catalog_ids: Vec<u64>,

    /// Pattern hash table (open addressing, quadratic probing).
    /// Each slot holds [star_idx0, star_idx1, star_idx2, star_idx3].
    /// Empty slots are [0, 0, 0, 0]. Since patterns always have 4 *distinct*
    /// star indices, all-zero is an unambiguous empty marker.
    pub pattern_catalog: Vec<[u32; 4]>,

    /// Largest edge angle (radians) for each pattern slot.
    /// Used for fast FOV-consistency filtering during solve.
    pub pattern_largest_edge: Vec<f32>,

    /// Lower 16 bits of the pattern key hash for each slot.
    /// Used as a fast pre-filter before full edge-ratio comparison.
    pub pattern_key_hashes: Vec<u16>,

    /// Database generation parameters.
    pub props: DatabaseProperties,
}

// ── Configuration for database generation ───────────────────────────────────

/// Parameters controlling database generation.
pub struct GenerateDatabaseConfig {
    /// Maximum FOV in degrees.
    pub max_fov_deg: f32,
    /// Minimum FOV in degrees. If None, equals max_fov (single-scale database).
    pub min_fov_deg: Option<f32>,
    /// Faintest star magnitude to include. None = auto-compute from star density.
    pub star_max_magnitude: Option<f32>,
    /// Maximum edge-ratio error for pattern matching.
    /// Determines bin count: bins = round(0.25 / pattern_max_error).
    /// Default 0.001 → 250 bins.
    pub pattern_max_error: f32,
    /// Oversampling factor for lattice field distribution. Default 100.
    pub lattice_field_oversampling: u32,
    /// Patterns to generate per lattice field. Default 50.
    pub patterns_per_lattice_field: u32,
    /// Max catalog stars per FOV for verification. Default 150.
    pub verification_stars_per_fov: u32,
    /// Multiscale FOV step ratio. Default 1.5.
    pub multiscale_step: f32,
    /// Year for proper motion propagation. None = don't propagate.
    pub epoch_proper_motion_year: Option<f64>,
    /// HEALPix nside for the spatial star catalog index. Default 16.
    pub catalog_nside: u32,
}

impl Default for GenerateDatabaseConfig {
    fn default() -> Self {
        Self {
            max_fov_deg: 30.0,
            min_fov_deg: None,
            star_max_magnitude: None,
            pattern_max_error: 0.001,
            lattice_field_oversampling: 100,
            patterns_per_lattice_field: 50,
            verification_stars_per_fov: 150,
            multiscale_step: 1.5,
            epoch_proper_motion_year: Some(2025.0),
            catalog_nside: 16,
        }
    }
}

// ── Configuration for plate solving ─────────────────────────────────────────

/// Parameters controlling the plate-solve attempt.
pub struct SolveConfig {
    /// Estimated field of view in radians.
    pub fov_estimate_rad: f32,
    /// Maximum acceptable FOV error in radians. None = no FOV filtering.
    pub fov_max_error_rad: Option<f32>,
    /// Maximum match distance as a fraction of the FOV. Default 0.01.
    pub match_radius: f32,
    /// False-positive probability threshold. Default 1e-5.
    pub match_threshold: f64,
    /// Timeout in milliseconds. None = no timeout. Default 5000.
    pub solve_timeout_ms: Option<u64>,
    /// Maximum edge-ratio error for matching. None = use database value.
    pub match_max_error: Option<f32>,
}

impl SolveConfig {
    /// Create a solve configuration with the given FOV estimate (radians).
    pub fn new(fov_estimate_rad: f32) -> Self {
        Self {
            fov_estimate_rad,
            fov_max_error_rad: None,
            match_radius: 0.01,
            match_threshold: 1e-5,
            solve_timeout_ms: Some(5000),
            match_max_error: None,
        }
    }
}

// ── Solve result ────────────────────────────────────────────────────────────

/// Result of a plate-solve attempt.
#[derive(Debug, Clone)]
pub struct SolveResult {
    /// Quaternion rotating ICRS vectors to camera-frame vectors.
    /// Camera frame: +X right, +Y down, +Z boresight.
    /// Usage: `camera_vec = quaternion * icrs_vec`
    pub quaternion: Option<Quaternion>,
    /// Estimated field of view (radians).
    pub fov_rad: Option<f32>,
    /// Number of matched stars in the verification step.
    pub num_matches: Option<u32>,
    /// RMS angular residual (radians) of matched stars.
    pub rmse_rad: Option<f32>,
    /// 90th-percentile angular residual (radians).
    pub p90e_rad: Option<f32>,
    /// Maximum angular residual (radians).
    pub max_err_rad: Option<f32>,
    /// False-positive probability (lower is better).
    pub prob: Option<f64>,
    /// Wall-clock time spent solving, in milliseconds.
    pub solve_time_ms: f32,
    /// Outcome status.
    pub status: SolveStatus,
    /// Catalog IDs of matched stars (only populated on success).
    pub matched_catalog_ids: Vec<u64>,
    /// Indices into the input centroid slice for each match.
    pub matched_centroid_indices: Vec<usize>,
}

impl SolveResult {
    /// Create a failure result with the given status and elapsed time.
    pub(crate) fn failure(status: SolveStatus, solve_time_ms: f32) -> Self {
        Self {
            quaternion: None,
            fov_rad: None,
            num_matches: None,
            rmse_rad: None,
            p90e_rad: None,
            max_err_rad: None,
            prob: None,
            solve_time_ms,
            status,
            matched_catalog_ids: Vec::new(),
            matched_centroid_indices: Vec::new(),
        }
    }
}