gem-hunter

A blind astrometry plate solver written in Rust, inspired by astrometry.net.

Given an astronomical image (or a list of detected source positions), zodiacal determines the sky coordinates — pointing, rotation, and plate scale — using only the pattern of stars in the image. No prior knowledge of the field is required.

Zodiacal implements the geometric hashing approach described in Lang et al. (2010):

1. Source extraction — detect bright point sources in the image
2. Quad formation — form 4-star asterisms from the brightest sources and compute scale-invariant hash codes
3. Code matching — search prebuilt index files for matching asterisms using a kd-tree
4. Hypothesis fitting — compute a TAN-WCS (tangent-plane) projection from matched star correspondences
5. Bayesian verification — evaluate the alignment against known catalog positions using a log-odds decision framework

The solver uses multi-scale indexes with per-index scale filtering to efficiently search across different angular scales.

Tested on 1,000 simulated fields (9568x6380 px, plate scale 0.1296"/px, FOV ~20'x13'):

Best configuration solve time distribution (985 solved):

• 48% solved in under 1 second
• 83% solved in under 5 seconds
• 93% solved in under 10 seconds
• 99% solved in under 30 seconds
• 60-second timeout for all cases

The 15 remaining failures are quad-matching timeouts where the solver exhausts the search budget without finding the correct asterism match. Analysis shows these are a mix of sparse fields (few detected sources) and fields where the correct quad is buried among too many candidates.

For comparison, the reference astrometry.net system reports 99.9% on SDSS fields (2048x2048 px, 0.396"/px, wider FOV ~13.5').

The 1000-field JSON corpus used for the table above lives in a separate data repo: OrbitalCommons/zodiacal-test-cases. Each NNNN.json file is a zodiacal::extraction::SourcesJson with truth RA/Dec, plate scale, and an extracted source list — drop-in input to zodiacal batch-solve and the zodiacal-tools bench-bundle harness.

Two sets are available:

• set1-legacy/ — the original meter-sim sensor_shootout corpus (the historical 96.7% baseline). Imported from this repo's old in-tree test_cases/.
• set2-dr3-mag19/ — regenerated against Gaia DR3 mag-19 via focalplane's motion_simulator.

Fetch a sibling checkout with the helper script:

scripts/fetch_test_cases.sh           # clones/updates ../zodiacal-test-cases
scripts/fetch_test_cases.sh --ssh     # use git@github.com remote
scripts/fetch_test_cases.sh --dest /some/path

The bench harness (and the analyze_failures.py / plot_failures.py scripts) default to ../zodiacal-test-cases/set1-legacy, so once the sibling checkout exists no further config is needed:

cargo run -p zodiacal-tools --release -- bench-bundle \
    --bundle-path <bundle.zdcl> \
    --test-cases-dir ../zodiacal-test-cases/set1-legacy

bench-bundle --trace-out DIR emits a <case>.trace.json sidecar per successful solve, containing the full WCS, the matched quad's four (field, catalog) pairs, and every verification matched-pair the solver recorded. Useful when you want to see which detection became A/B/C/D and which top-50 stars verified vs missed.

cargo run -p zodiacal-tools --release -- bench-bundle \
    --bundle-path <bundle.zdcl> \
    --test-cases-dir ../zodiacal-test-cases/set2-dr3-mag19 \
    --trace-out /tmp/traces

scripts/render_solver_trace.py overlays a trace on its FITS plate using the manim default palette (A=red, B=green, C=cyan, D=orange), with clipped backbone/reference lines and a dashed AB-diameter validity circle. Top-50 verification hits draw as open lime circles, misses as crimson ×s, and rank > 50 detections as faint cyan ×s underneath. Requires numpy, matplotlib, and astropy.

python scripts/render_solver_trace.py \
    --case 0042 \
    --fits  /path/to/frame_0042.fits \
    --src-json  ../zodiacal-test-cases/set2-dr3-mag19/0042.json \
    --trace /tmp/traces/0042.trace.json \
    --out   /tmp/0042.png

Zodiacal requires prebuilt index files (.zdcl) containing star positions and quad hash codes. These are built from a starfield binary catalog.

Star catalogs are built using the gaia_filter tool from the starfield crate, which filters cached Gaia DR3 source files by magnitude:

# Build a magnitude-19 catalog (~565M stars, ~17 GB)
cd /path/to/starfield
cargo build --release --example gaia_filter
./target/release/examples/gaia_filter --magnitude 19.0 --output gaia_mag19.bin --threads 8 all

The all argument processes all cached Gaia source files. Brighter magnitude limits produce smaller, faster-to-load catalogs:

The recommended approach builds a series of narrow-band indexes, each covering a different angular scale range:

zodiacal build-index-series \
  -c gaia_mag19.bin \
  -o my_index \
  --scale-lower 10.0 \
  --scale-upper 600.0 \
  --max-stars-per-cell 30 \
  --passes 16 \
  --max-reuse 8 \
  --max-depth 8

This produces index files my_index_00.zdcl through my_index_11.zdcl, with 12 bands at a sqrt(2) scale factor. Each band covers a narrow range of quad angular sizes (e.g., 80-113", 113-160"), keeping the code trees small to reduce false-positive matches.

Key parameters:

• --scale-lower/--scale-upper: angular scale range in arcseconds
• --scale-factor: ratio between adjacent bands (default sqrt(2) ~1.414; use 3.0 for fewer, wider bands)
• --max-stars-per-cell: brightest stars per HEALPix cell for uniformization (default 30)
• --max-depth: maximum HEALPix depth (default 8 = 786,432 cells; limits memory usage)
• --passes: number of quad-building passes per cell
• --max-reuse: maximum times a star can be reused across quads

For quick testing or narrow-scale applications:

zodiacal build-index \
  -c gaia_mag16.bin \
  -o test_index.zdcl \
  --scale-lower 30.0 \
  --scale-upper 600.0

zodiacal solve image.png \
  -i my_index_00.zdcl \
  -i my_index_01.zdcl \
  -i my_index_02.zdcl \
  --timeout 60

Process a directory of images or JSON source files:

zodiacal batch-solve images/ \
  -i my_index_00.zdcl -i my_index_01.zdcl -i my_index_02.zdcl \
  --pattern "*.json" \
  --timeout 60 \
  --failures-dir /tmp/failures

The --failures-dir option writes diagnostic JSON for each unsolved case, useful for analyzing failure modes.

The solver accepts JSON source lists directly (skips extraction):

zodiacal solve sources.json \
  -i my_index_00.zdcl -i my_index_01.zdcl

Extract sources from an image and write them as JSON for later solving:

zodiacal extract image.png -o sources.json

Zodiacal uses a JSON format for exchanging detected source lists between tools. The extract command produces this format, and the solver can consume it.

{
  "image_width": 9568,
  "image_height": 6380,
  "ra_deg": 265.47,
  "dec_deg": 44.31,
  "plate_scale_arcsec": 0.1296,
  "sources": [
    { "x": 4821.3, "y": 3190.7, "flux": 54210.0 },
    { "x": 1023.5, "y": 892.1, "flux": 38450.0 }
  ]
}

The ra_deg, dec_deg, and plate_scale_arcsec fields are optional hints. When present they can speed up solving by constraining the search space. When absent they are omitted from the JSON entirely.

The CLI is a thin wrapper over the zodiacal crate. Embedding the solver into your own application gives you access to several features the CLI doesn't expose, plus three deployment-mode building blocks (see below).

use zodiacal::index::Index;
use zodiacal::solver::{solve, SolverConfig};
use zodiacal::extraction::DetectedSource;

let index = Index::load("my_index_00.zdcl".as_ref())?;
let sources: Vec<DetectedSource> = /* from extract_sources() or your own pipeline */;
let (solution, stats) = solve(&sources, &[&index], (image_w, image_h), &SolverConfig::default());

zodiacal builds for wasm32-unknown-unknown with default features off:

zodiacal = { version = "0.5", default-features = false }

Solving, verification, fitting and in-memory index building work in the browser, and solver timeouts use the browser clock. Anything that touches the filesystem does not: Index::load, LiveIndex, memory-mapped index sources, index and sidecar writers, and starfield's datastore cache are native-only. In the browser, fetch index data in the host and construct an Index in memory (its fields are public, or use index::builder::build_index on a star list). rayon runs single-threaded. CI checks the wasm32 build and runs solver tests under Node (tests/wasm.rs).

After the blind solve produces a TAN+SIP solution, the refinement module re-fits the WCS using each matched catalog star's apparent direction at the observation time — applying proper motion, parallax, light-time, and stellar aberration via the starfield apparent-place pipeline.

use zodiacal::refinement::{refine_solution, ObservationContext, ObserverState, RefinementConfig};

let obs = ObservationContext {
    time: starfield::time::Timescale::default().tt_jd(jd_tt, None),
    observer: ObserverState::Barycentric {
        position_au: spacecraft_pos,
        velocity_au_per_day: spacecraft_vel,
    },
};
let refined = refine_solution(&tweaked_wcs, &sources, &index, &gaia_catalog, &obs, &RefinementConfig::default())?;
println!("residual RMS: {:.2} mas", refined.residual_rms_mas);

The full Gaia astrometry needed for refinement (RA/Dec/PM/parallax/RV + per-component sigmas) lives in a RefinementCatalog. For real datasets you'd populate it from the Gaia sidecar (my_index.zdcl.gaia) — a flat sorted-by-source_id binary file colocated with the index, accessed via mmap + galloping search:

use zodiacal::refinement::RefinementCatalog;

let catalog_ids: Vec<u64> = matched_sources.iter().map(|s| s.catalog_id).collect();
let gaia = RefinementCatalog::load_sidecar_filtered(
    "my_index_00.zdcl.gaia".as_ref(),
    &catalog_ids,
)?;

The library exposes three building blocks that compose into different operational profiles. The summaries below get you started; for end-to-end code, memory budgets, performance numbers, sharding, and pitfalls per mode, see the dedicated docs:

• docs/deployment/server.md — backend / batch service
• docs/deployment/ground.md — ground telescope tracking the zenith
• docs/deployment/space.md — spacecraft star tracker
• docs/deployment/README.md — overview + cross-cutting performance table

Load every index once, share across many concurrent solves. The standard Index::load + solve API works as-is. Wrap in Arc<Index> to share across threads.

use std::sync::Arc;
let index = Arc::new(Index::load("my_index_00.zdcl".as_ref())?);
// dispatch solves across rayon / tokio / thread pool

A ground telescope's visible sky changes as Earth rotates. LiveIndex tracks which HEALPix cells of the index need to be resident, GroundStation provides the current zenith from (lat, lon, time), and RealtimeSolver ticks them in sync:

use zodiacal::index::{LiveIndex, ZdclFile};
use zodiacal::pointing::GroundStation;
use zodiacal::realtime::{RealtimeSolver, RefreshPolicy};
use std::time::Duration;

let source = ZdclFile::open("my_index_v3.zdcl".as_ref())?;
let pointing = GroundStation::from_degrees(34.0, -118.0, /* min altitude */ 30.0);
let mut rt = RealtimeSolver::new(source, pointing)
    .with_refresh_policy(RefreshPolicy::OnPointingDelta {
        angular_threshold_rad: 0.1,
        max_age: Duration::from_secs(60),
    });
// later, per frame:
let out = rt.solve(&detected_sources, image_size, &now)?;

Same orchestrator, with a SpacecraftBoresight driving the loaded region from an attitude estimate + ephemeris. Bring your own EphemerisSource (anise, SPICE, TLE via starfield::sgp4lib, custom Kalman) and AttitudeSource:

use zodiacal::pointing::{SpacecraftBoresight, EphemerisSource, AttitudeSource};

let pointing = SpacecraftBoresight::new(my_ephemeris, my_attitude_estimator,
    /* detector half-angle */ 0.05_f64.to_radians(),
    /* fov padding */ 0.05_f64.to_radians());
let mut rt = RealtimeSolver::new(zdcl_file_v3, pointing);

Spacecraft mode also feeds observer_state(t) into refinement automatically — you get parallax + aberration corrections for free using the same ephemeris.

For server-mode callers with prior knowledge of where they're pointing, Index::load_in_region reads the same v2 file format but keeps only stars inside a SkyRegion. Same disk I/O as a full load, dramatically lower resident memory.

use zodiacal::solver::SkyRegion;
use starfield::Equatorial;
let region = SkyRegion::from_degrees(Equatorial::new(0.5, 0.3), 5.0);
let small_index = Index::load_in_region_padded(path, &region, /* pad rad */ 0.01)?;

.zdcl index files are versioned. Both the original streaming format and the new HEALPix-grouped layout are accepted by all current readers.

ZdclFile::open accepts v1/v2/v3 transparently; older files appear as a single virtual cell so the IndexSource API works uniformly.

The blind solve pipeline:

Image → Source Extraction → Quad Formation → Code Matching → WCS Fitting → Verification
                                                  ↑
                                          Index Files (.zdcl)
                                                  ↑
                                     Star Catalog → HEALPix Uniformization → Quad Building

Optional post-solve refinement chain:

Tweaked WCS + Matched Sources + Gaia sidecar + Observation Context
                                  ↓
                         Apparent place per matched star
                         (PM + parallax + light-time + aberration)
                                  ↓
                         Weighted re-fit of TAN parameters
                                  ↓
                          Refined Solution (~10 mas RMS)

For realtime modes, LiveIndex sits between the file and the solver, holding only the HEALPix cells currently relevant. RealtimeSolver coordinates pointing-source updates with LiveIndex membership changes and the cached snapshot used for solving.

Key design decisions:

• HEALPix uniformization at index build, HEALPix grouping in v3 file format: stars are distributed across sky cells to ensure uniform coverage and to enable cell-targeted reads.
• Multi-scale indexes: narrow angular-scale bands keep kd-trees small, reducing false-positive matches.
• Dual-parity matching: field codes are tried in both orientations (original + x/y swapped) to handle image flips.
• Pre-fit filters: SolverConfig::scale_range and SolverConfig::within reject candidates before the LSQ fit, eliminating ~10× pessimization on wide hints.
• KdForest for live indexes: per-cell sub-trees so cell add/drop is O(1) — no full rebuild on every realtime membership change.
• Generation-cached snapshots: RealtimeSolver caches the flat Index view by LiveIndex::build_generation, so steady-state solve cost is a pointer-deref.

• starfield — star catalogs, coordinate systems, time scales, apparent-place pipeline
• cdshealpix — HEALPix tessellation for index uniformization and v3 cell layout
• memmap2 — mmap-backed sidecar and v3 index readers
• nalgebra — quaternion + vector math (used by refinement and pointing-source)
• ndarray (optional, image-processing feature) — N-dimensional arrays
• clap (optional, cli feature) — CLI argument parsing
• image (optional, cli feature) — image loading
• rayon — parallelism
• serde / serde_json — serialization

See LICENSE for details.