arcsec_core/pipeline/

solver.rs

//! Full plate-solving pipeline.

use core::f64::consts::PI;
use std::path::PathBuf;

use crate::catalog::read_catalog_stars;
use crate::catalog::{CatalogLayout, CatalogStar};
use crate::detection::get_background;
use crate::detection::stars::find_stars_with_background;
use crate::error::{ArcsecError, Result};
use crate::math::coords::{ang_sep, equatorial_standard, standard_equatorial};
use crate::math::lsq::{fit_affine, solve_plate_constants};
use crate::quads::{
    TETRA_TOL_FACTOR, bijective_filter, build_quads, build_quads_presorted, build_triangles,
    extract_star_pairs, extract_triangle_pairs, filter_by_scale, filter_triangles_by_scale,
    find_matches_sorted, find_triangle_matches, vote_filter,
};
use crate::types::{MatchedStar, PairedPositions, PlateConstants, Star, StarList, WcsSolution};
use crate::wcs::output::derive_wcs;

use super::spiral::SpiralSearch;

/// Which pattern-matching algorithm to use in the catalog spiral loop.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum SolveMethod {
    /// ASTAP-style 5-ratio quad matching with `vote_filter` (default).
    #[default]
    Quads,
    /// TETRA 2-ratio triangle matching with bijective filter.
    Tetra,
}

/// How much sky the spiral search reads around each position (ASTAP's `-speed`).
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum SearchSpeed {
    /// Size the catalogue window by the image's star count: twice the field for an
    /// image with fewer than 35 stars, falling to the field itself above 140.
    #[default]
    Auto,
    /// Always read a window twice the field, so neighbouring spiral positions
    /// overlap and a field that straddles two of them is still seen whole. Four
    /// times the catalogue stars per position, so slower; it helps images with
    /// many stars that the auto window only just misses.
    Slow,
}

/// Parameters for [`solve_image`].
#[derive(Debug, Clone)]
pub struct SolveParams {
    /// Approximate RA of image centre (radians, hint only).
    pub ra_hint: f64,
    /// Approximate DEC of image centre (radians, hint only).
    pub dec_hint: f64,
    /// Image field of view (square side, radians). Used as the spiral step size.
    pub fov: f64,
    /// Maximum search radius from the hint position (radians).
    pub search_radius: f64,
    /// Quad ratio matching tolerance (ASTAP default ≈ 0.007).
    pub quad_tolerance: f64,
    /// Minimum HFD for valid stars (pixels).
    pub hfd_min: f64,
    /// Maximum number of image stars to detect.
    pub max_stars: usize,
    /// Path to the catalog directory.
    pub db_path: PathBuf,
    /// Catalog name prefix (e.g. `"d20"`, `"d80"`).
    pub db_name: String,
    /// Pixel binning factor applied before solving (1 = none, 2 = 2×2, ...).
    /// WCS output is scaled back to original image pixel coordinates.
    pub binning: usize,
    /// Pattern-matching algorithm for the catalog spiral.
    pub method: SolveMethod,
    /// Catalogue window per spiral position.
    pub speed: SearchSpeed,
    /// Worker threads for the spiral search. 0 = one per available core.
    ///
    /// Spiral positions are independent, so they are evaluated a batch at a time
    /// across this many threads. Results are identical to the serial search: within
    /// a batch the lowest spiral index still wins, so the first position that
    /// verifies is the one returned, exactly as before.
    pub threads: usize,
}

/// Iterative sigma-clipping: fit plate constants, reject pairs with large residuals,
/// re-fit until stable or fewer than `min_count` pairs remain.
///
/// First pass uses a 10-pixel absolute threshold (in catalog arcsec) to cut the
/// large residuals of false-positive triangle matches. Subsequent passes apply
/// `sigma × rms` clipping until the set is stable.
fn sigma_clip_pairs(
    mut img_pos: Vec<(f64, f64)>,
    mut cat_pos: Vec<(f64, f64)>,
    sigma: f64,
    min_count: usize,
) -> PairedPositions {
    let mut first_pass = true;
    for _ in 0..10 {
        if img_pos.len() < min_count.max(3) {
            break;
        }
        // Unchecked: the first fit is made on the contaminated set, and gross
        // outliers can skew it past solve_plate_constants' scale check even though
        // clipping them is exactly what would fix it.
        let Ok(plate) = fit_affine(&img_pos, &cat_pos) else {
            break;
        };
        let residuals: Vec<f64> = img_pos
            .iter()
            .zip(cat_pos.iter())
            .map(|(&(xi, yi), &(xc, yc))| {
                let xp = plate.a * xi + plate.b * yi + plate.c;
                let yp = plate.d * xi + plate.e * yi + plate.f;
                ((xp - xc).powi(2) + (yp - yc).powi(2)).sqrt()
            })
            .collect();
        let rms = (residuals.iter().map(|r| r * r).sum::<f64>() / residuals.len() as f64).sqrt();
        let threshold = if first_pass {
            first_pass = false;
            // 10 px in catalog-arcsec: generous cut for large FP residuals on first pass.
            let cdelt = (plate.a.powi(2) + plate.d.powi(2)).sqrt();
            // Gross outliers drag a least-squares fit towards themselves and inflate
            // every residual, so a fixed cut can keep them. The median residual is
            // not moved by a minority of outliers: allow 3 sigma of it (1.4826 x the
            // median absolute residual estimates sigma) when that is larger.
            let mut sorted = residuals.clone();
            sorted.sort_unstable_by(f64::total_cmp);
            let median = sorted[sorted.len() / 2];
            (10.0 * cdelt).max(10.0).max(3.0 * 1.4826 * median)
        } else {
            sigma * rms
        };
        let before = img_pos.len();
        let mut new_img = Vec::with_capacity(before);
        let mut new_cat = Vec::with_capacity(before);
        for ((&ip, &cp), &r) in img_pos.iter().zip(cat_pos.iter()).zip(residuals.iter()) {
            if r <= threshold {
                new_img.push(ip);
                new_cat.push(cp);
            }
        }
        if new_img.len() == before {
            break; // stable — no more outliers
        }
        img_pos = new_img;
        cat_pos = new_cat;
    }
    (img_pos, cat_pos)
}

/// Minimum number of individually matched stars required to believe a solution.
///
/// Correct solves typically match 200-375 stars, so this is deliberately loose;
/// its job is to reject the handful-of-coincidences case. Together with
/// `MIN_VERIFY_SPREAD` it separates two otherwise identical-looking results: M31 at
/// 2 degrees (22 stars, spread 0.221, rms 0.65", rotation wrong by 1.56 degrees)
/// from the Dec -88 field (46 stars, spread 0.207, rms 0.66", correct to 2.3").
const MIN_VERIFIED_STARS: usize = 30;
/// Match radii (pixels) used by successive verification passes, coarse to fine.
const VERIFY_RADII: [f64; 3] = [6.0, 3.0, 2.0];
/// Minimum spread of the matched stars, as a fraction of the image half-diagonal.
///
/// A count threshold alone is not enough: matches clustered in one part of the
/// frame (the core of a bright galaxy, say) pin the position but leave rotation
/// and scale essentially free. M31 at 2 degrees passed with 22 matched stars and
/// a 1.56-degree rotation error, which is 154" at the field corners.
const MIN_VERIFY_SPREAD: f64 = 0.20;

/// A verified plate: the re-fitted plate constants, the per-star RMS in arcsec,
/// and the star pairs the fit was made from.
struct Verified {
    plate: PlateConstants,
    rms: f64,
    /// Detected star positions, pixels of the solved (binned) image, 0-based.
    img_pos: Vec<(f64, f64)>,
    /// The catalogue star each was paired with, in standard coordinates (arcsec)
    /// about the plane the plate maps into.
    cat_pos: Vec<(f64, f64)>,
}

impl Verified {
    /// Number of individually matched stars.
    fn n(&self) -> usize {
        self.img_pos.len()
    }
}

/// Project the catalogue onto the image with a candidate plate solution, match
/// individual stars, and re-fit on those matches.
///
/// The quad matcher only ever produces quad *centroids*, so the plate fit is built
/// from a handful of averaged positions and nothing ever checks that the individual
/// stars agree. This does that check: invert the plate to map every catalogue star
/// into pixel space, pair each with the nearest detected star, re-fit on the pairs,
/// and repeat with a shrinking radius.
///
/// Returns the refined plate with its matched pairs, or `None` if the plate is
/// degenerate, too few stars agree, or the matches are too clustered.
fn verify_and_refit(
    img_stars: &StarList,
    cat_stars: &StarList,
    plate: &PlateConstants,
    img_w: usize,
    img_h: usize,
) -> Option<Verified> {
    if img_stars.is_empty() || cat_stars.is_empty() {
        return None;
    }

    // Uniform grid over the detected stars for nearest-neighbour lookup.
    let (mut min_x, mut min_y) = (f64::INFINITY, f64::INFINITY);
    let (mut max_x, mut max_y) = (f64::NEG_INFINITY, f64::NEG_INFINITY);
    for st in &img_stars.0 {
        min_x = min_x.min(st.x);
        max_x = max_x.max(st.x);
        min_y = min_y.min(st.y);
        max_y = max_y.max(st.y);
    }
    if !(min_x.is_finite() && min_y.is_finite() && max_x > min_x && max_y > min_y) {
        return None;
    }
    let cell = VERIFY_RADII[0].max(1.0);
    let nx = (((max_x - min_x) / cell).ceil() as usize + 1).max(1);
    let ny = (((max_y - min_y) / cell).ceil() as usize + 1).max(1);
    let mut grid: Vec<Vec<u32>> = vec![Vec::new(); nx * ny];
    for (i, st) in img_stars.0.iter().enumerate() {
        let gx = ((st.x - min_x) / cell) as usize;
        let gy = ((st.y - min_y) / cell) as usize;
        grid[gy.min(ny - 1) * nx + gx.min(nx - 1)].push(i as u32);
    }

    let mut current = plate.clone();
    // The last pass that fitted, with the spread of its matches.
    let mut best: Option<(Verified, f64)> = None;

    for &radius in &VERIFY_RADII {
        let det = current.a * current.e - current.b * current.d;
        if det.abs() < 1e-12 {
            return None;
        }
        let r2 = radius * radius;

        let mut img_pos: Vec<(f64, f64)> = Vec::new();
        let mut cat_pos: Vec<(f64, f64)> = Vec::new();
        let mut used = vec![false; img_stars.len()];

        for cs in &cat_stars.0 {
            // Invert  xi = a*x + b*y + c ;  eta = d*x + e*y + f
            let dx = cs.x - current.c;
            let dy = cs.y - current.f;
            let px = (current.e * dx - current.b * dy) / det;
            let py = (-current.d * dx + current.a * dy) / det;
            if px < min_x - radius
                || px > max_x + radius
                || py < min_y - radius
                || py > max_y + radius
            {
                continue;
            }

            let gx = (((px - min_x) / cell) as isize).clamp(0, nx as isize - 1);
            let gy = (((py - min_y) / cell) as isize).clamp(0, ny as isize - 1);
            let mut best_i: Option<usize> = None;
            let mut best_d2 = r2;
            for oy in -1isize..=1 {
                for ox in -1isize..=1 {
                    let cx = gx + ox;
                    let cy = gy + oy;
                    if cx < 0 || cy < 0 || cx >= nx as isize || cy >= ny as isize {
                        continue;
                    }
                    for &i in &grid[cy as usize * nx + cx as usize] {
                        let i = i as usize;
                        if used[i] {
                            continue;
                        }
                        let st = &img_stars.0[i];
                        let d2 = (st.x - px) * (st.x - px) + (st.y - py) * (st.y - py);
                        if d2 < best_d2 {
                            best_d2 = d2;
                            best_i = Some(i);
                        }
                    }
                }
            }
            if let Some(i) = best_i {
                used[i] = true; // one-to-one: a detected star backs at most one catalogue star
                img_pos.push((img_stars.0[i].x, img_stars.0[i].y));
                cat_pos.push((cs.x, cs.y));
            }
        }

        if img_pos.len() < 4 {
            break;
        }
        let Ok(refined) = solve_plate_constants(&img_pos, &cat_pos) else {
            break;
        };
        let mut sq = 0.0;
        for (&(xi, yi), &(xc, yc)) in img_pos.iter().zip(cat_pos.iter()) {
            let xp = refined.a * xi + refined.b * yi + refined.c;
            let yp = refined.d * xi + refined.e * yi + refined.f;
            sq += (xp - xc).powi(2) + (yp - yc).powi(2);
        }
        let rms = (sq / img_pos.len() as f64).sqrt();
        // Spread of the matched stars about their own centroid, as a fraction of the
        // image half-diagonal. Matches clustered in one corner leave rotation free.
        let n = img_pos.len() as f64;
        let mx = img_pos.iter().map(|p| p.0).sum::<f64>() / n;
        let my = img_pos.iter().map(|p| p.1).sum::<f64>() / n;
        let var = img_pos
            .iter()
            .map(|&(x, y)| (x - mx) * (x - mx) + (y - my) * (y - my))
            .sum::<f64>()
            / n;
        let half_diag = 0.5 * ((img_w * img_w + img_h * img_h) as f64).sqrt();
        let spread = var.sqrt() / half_diag;
        log::debug!(
            "verify: {} stars, spread {:.3}, rms {:.2}\"",
            img_pos.len(),
            spread,
            rms
        );

        current = refined.clone();
        best = Some((
            Verified {
                plate: refined,
                rms,
                img_pos,
                cat_pos,
            },
            spread,
        ));
    }

    best.filter(|(v, spread)| v.n() >= MIN_VERIFIED_STARS && *spread >= MIN_VERIFY_SPREAD)
        .map(|(v, _)| v)
}

/// Everything a spiral position needs that does not change between positions.
struct SpiralCtx<'a> {
    params: &'a SolveParams,
    img: &'a crate::types::ImageBuffer,
    stars: &'a StarList,
    img_quads: &'a crate::types::QuadList,
    img_tris: &'a crate::quads::TriangleList,
    nrstars_image: usize,
    nrstars_required: usize,
    oversize: f64,
    min_quads: usize,
    step_size: f64,
}

/// A spiral position that produced a verified solution.
struct PositionOutcome {
    idx: usize,
    ra_db: f64,
    dec_db: f64,
    sep_deg: f64,
    verified: Verified,
    n_matched: usize,
    n_raw: usize,
    mag_limit: f64,
}

/// Result of trying one spiral position: the angular distance if the catalogue was
/// actually read there (for the ASTAP-style progress line), and the solution if one
/// verified.
struct PositionTry {
    sep_deg: Option<f64>,
    outcome: Option<PositionOutcome>,
}

impl PositionTry {
    const NONE: Self = Self {
        sep_deg: None,
        outcome: None,
    };
}

/// Evaluate a single spiral position. Pure with respect to `ctx`, so positions can
/// be run concurrently.
fn try_position(ctx: &SpiralCtx<'_>, idx: usize, sx: i32, sy: i32) -> PositionTry {
    let params = ctx.params;
    let step_size = ctx.step_size;

    let dec_db_raw = params.dec_hint + step_size * sy as f64;
    let (dec_db, flip) = if dec_db_raw > PI / 2.0 {
        (PI - dec_db_raw, PI)
    } else if dec_db_raw < -PI / 2.0 {
        (-PI - dec_db_raw, PI)
    } else {
        (dec_db_raw, 0.0)
    };

    let extra = if dec_db > 0.0 {
        step_size * 0.5
    } else {
        -step_size * 0.5
    };
    let ra_offset = step_size * sx as f64 / (dec_db - extra).cos();
    if ra_offset > PI / 2.0 + step_size * 0.5 || ra_offset < -PI / 2.0 {
        return PositionTry::NONE;
    }

    let ra_db = (flip + params.ra_hint + ra_offset).rem_euclid(2.0 * PI);
    let sep = ang_sep(ra_db, dec_db, params.ra_hint, params.dec_hint);
    if sep > params.search_radius + step_size / 2.0 {
        return PositionTry::NONE;
    }

    // Any read failure (a missing tile included) counts as "nothing catalogued
    // here"; `solve_image` has already checked that the database exists at all.
    let cat_raw = match read_catalog_stars(
        &params.db_path,
        &params.db_name,
        ra_db,
        dec_db,
        params.fov * ctx.oversize,
        ctx.nrstars_required,
    ) {
        Ok(v) if !v.is_empty() => v,
        Ok(_) | Err(_) => return PositionTry::NONE,
    };

    let sep_deg = sep.to_degrees();
    let mag_limit = cat_raw
        .iter()
        .map(|s| s.mag)
        .fold(f64::NEG_INFINITY, f64::max);
    log::info!(
        "Search {}, [{},{}], position: {}  Down to magn {:.1}  {} database stars  {} database quads to compare.",
        idx,
        sx,
        sy,
        format_radec(ra_db, dec_db),
        mag_limit,
        cat_raw.len(),
        cat_raw.len(),
    );

    let mut cat_stars: Vec<Star> = cat_raw
        .iter()
        .map(|s| {
            let (x, y) = equatorial_standard(ra_db, dec_db, s.ra, s.dec, 1.0);
            Star {
                x,
                y,
                snr: 1.0,
                hfd: 2.0,
            }
        })
        .collect();
    cat_stars.sort_unstable_by(|a, b| a.x.total_cmp(&b.x));
    let cat_star_list = StarList(cat_stars);

    let failed = PositionTry {
        sep_deg: Some(sep_deg),
        outcome: None,
    };

    let (img_pos, cat_pos, n_matched, n_raw) = match params.method {
        SolveMethod::Quads => {
            let mut cat_quads = build_quads_presorted(&cat_star_list, ctx.nrstars_image);
            if cat_quads.is_empty() {
                return failed;
            }
            crate::quads::r#match::sort_catalog_quads(&mut cat_quads);
            let raw = find_matches_sorted(ctx.img_quads, &cat_quads, params.quad_tolerance);
            let n_raw = raw.len();
            log::info!("Found {n_raw} references");
            let mut filtered = vote_filter(ctx.img_quads, &cat_quads, &raw, params.quad_tolerance);
            if filtered.len() < ctx.min_quads {
                let (by_scale, _) = filter_by_scale(&raw, params.quad_tolerance);
                if by_scale.len() > filtered.len() {
                    filtered = by_scale;
                }
            }
            if filtered.len() < ctx.min_quads {
                return failed;
            }
            let (ip, cp) = extract_star_pairs(ctx.img_quads, &cat_quads, &filtered);
            (ip, cp, filtered.len(), n_raw)
        }
        SolveMethod::Tetra => {
            let cat_tris = build_triangles(&cat_star_list);
            if cat_tris.is_empty() {
                return failed;
            }
            let tol = params.quad_tolerance * TETRA_TOL_FACTOR;
            let raw = find_triangle_matches(ctx.img_tris, &cat_tris, tol);
            let n_raw = raw.len();
            log::info!("Found {n_raw} triangle references");
            let biject = bijective_filter(&raw, ctx.img_tris, &cat_tris);
            let (filtered, _) = filter_triangles_by_scale(&biject, params.quad_tolerance);
            if filtered.len() < ctx.min_quads {
                return failed;
            }
            let (ip, cp) = extract_triangle_pairs(ctx.img_tris, &cat_tris, &filtered);
            let (ip, cp) = sigma_clip_pairs(ip, cp, 3.0, ctx.min_quads);
            if ip.len() < ctx.min_quads {
                return failed;
            }
            let n_clean = ip.len();
            (ip, cp, n_clean, n_raw)
        }
    };

    let Ok(plate) = solve_plate_constants(&img_pos, &cat_pos) else {
        return failed;
    };

    let Some(verified) = verify_and_refit(
        ctx.stars,
        &cat_star_list,
        &plate,
        ctx.img.width,
        ctx.img.height,
    ) else {
        log::info!("Verification failed at this position; continuing search.");
        return failed;
    };
    log::info!(
        "Verified {} stars against the catalogue, residual {:.2}\"",
        verified.n(),
        verified.rms
    );

    let (verified, ra_db, dec_db) = recentre(ctx, &cat_raw, verified, ra_db, dec_db);

    PositionTry {
        sep_deg: Some(sep_deg),
        outcome: Some(PositionOutcome {
            idx,
            ra_db,
            dec_db,
            sep_deg,
            verified,
            n_matched,
            n_raw,
            mag_limit,
        }),
    }
}

/// Refit a verified plate in the tangent plane at the image centre.
///
/// The plate constants are a linear map from pixels to the tangent plane at the
/// spiral position the catalogue was projected about. Pixels map linearly onto a
/// tangent plane only at the optical axis, so away from it the fit absorbs the
/// projection's curvature as a rotation and shear, which grow with the distance
/// from the field and with declination. Star-level RMS stays small, because the fit
/// is good *in that plane*, but the CD matrix derived from it is wrong at the image
/// centre: with the hint 0.3 fields off, most corpus solves were out by 5-1600" at
/// the corners.
///
/// So once a position verifies, move the tangent point to the image centre, pair
/// stars as the verified plate predicts them, fit those pairs in the new plane,
/// and verify again. Twice, since the centre moves slightly with
/// the new fit. If a pass fails to verify, the previous solution is kept: this can
/// only improve a solve, never lose one.
fn recentre(
    ctx: &SpiralCtx<'_>,
    cat_raw: &[CatalogStar],
    mut verified: Verified,
    mut ra_db: f64,
    mut dec_db: f64,
) -> (Verified, f64, f64) {
    let (w, h) = (ctx.img.width as f64, ctx.img.height as f64);
    let (cx, cy) = ((w - 1.0) * 0.5, (h - 1.0) * 0.5);
    let apply =
        |p: &PlateConstants, x: f64, y: f64| (p.a * x + p.b * y + p.c, p.d * x + p.e * y + p.f);

    for _ in 0..2 {
        let plate = &verified.plate;
        let (xs, ys) = apply(plate, cx, cy);
        // Already centred to well under a milliarcsecond: nothing to gain.
        if xs.hypot(ys) < 1e-3 {
            break;
        }
        let (ra0, dec0) = standard_equatorial(ra_db, dec_db, xs, ys, 1.0);

        // Starting plate for the new tangent plane. Mapping one tangent plane onto
        // another is far from linear over a wide field (a 10-degree field 3 degrees
        // off moves by ~65" under a straight-line fit), so rather than carry the
        // plate across, pair stars exactly as the verified plate predicts them and
        // fit those pairs against their positions in the new plane.
        let det = plate.a * plate.e - plate.b * plate.d;
        if det.abs() < 1e-12 {
            break;
        }
        let r2 = VERIFY_RADII[0] * VERIFY_RADII[0];
        let mut used = vec![false; ctx.stars.len()];
        let mut img_pos = Vec::new();
        let mut new_pos = Vec::new();
        let mut cat = Vec::with_capacity(cat_raw.len());
        for s in cat_raw {
            let (nx, ny) = equatorial_standard(ra0, dec0, s.ra, s.dec, 1.0);
            cat.push(Star {
                x: nx,
                y: ny,
                snr: 1.0,
                hfd: 2.0,
            });
            let (ox, oy) = equatorial_standard(ra_db, dec_db, s.ra, s.dec, 1.0);
            let (dx, dy) = (ox - plate.c, oy - plate.f);
            let px = (plate.e * dx - plate.b * dy) / det;
            let py = (-plate.d * dx + plate.a * dy) / det;
            let nearest = ctx
                .stars
                .0
                .iter()
                .enumerate()
                .filter(|&(i, _)| !used[i])
                .map(|(i, st)| (i, (st.x - px).powi(2) + (st.y - py).powi(2)))
                .filter(|&(_, d2)| d2 < r2)
                .min_by(|a, b| a.1.total_cmp(&b.1));
            if let Some((i, _)) = nearest {
                used[i] = true;
                img_pos.push((ctx.stars.0[i].x, ctx.stars.0[i].y));
                new_pos.push((nx, ny));
            }
        }
        let Ok(guess) = solve_plate_constants(&img_pos, &new_pos) else {
            break;
        };
        let cat = StarList(cat);
        let Some(v) = verify_and_refit(ctx.stars, &cat, &guess, ctx.img.width, ctx.img.height)
        else {
            log::info!("Re-centring on the image centre did not verify; keeping the fit.");
            break;
        };
        log::info!(
            "Re-centred on the image centre: verified {} stars, residual {:.2}\"",
            v.n(),
            v.rms
        );
        (verified, ra_db, dec_db) = (v, ra0, dec0);
    }
    (verified, ra_db, dec_db)
}

/// Solve the WCS for an image against an ASTAP star database.
///
/// Walks a square spiral out from the hint in steps of one field of view, and
/// returns the first position whose quad match survives star-by-star verification.
/// If `params.binning > 1`, `img` is taken to be the binned image and the returned
/// CRPIX/CD/CDELT are scaled back to the unbinned pixel grid.
///
/// All progress is emitted via the `log` crate at INFO level — callers install
/// whichever logger backend they need (file, stderr, both, or none).
///
/// # Errors
///
/// - [`ArcsecError::InvalidParameter`] if `fov` is not positive and finite, or
///   `search_radius` is negative or not finite.
/// - [`ArcsecError::CatalogNotFound`] if `db_path` holds no database called `db_name`.
/// - [`ArcsecError::InsufficientStars`] if fewer than 5 stars are detected.
/// - [`ArcsecError::InsufficientQuads`] if no spiral position yields a verified match.
pub fn solve_image(img: &crate::types::ImageBuffer, params: &SolveParams) -> Result<WcsSolution> {
    // The spiral steps by one FOV out to the search radius, so a zero, negative or
    // NaN FOV would make the step count infinite (and saturate to i32::MAX).
    if !(params.fov.is_finite() && params.fov > 0.0) {
        return Err(ArcsecError::InvalidParameter(format!(
            "field of view must be positive, got {} rad",
            params.fov
        )));
    }
    if !(params.search_radius.is_finite() && params.search_radius >= 0.0) {
        return Err(ArcsecError::InvalidParameter(format!(
            "search radius must be non-negative, got {} rad",
            params.search_radius
        )));
    }

    // Check the database up front. Every spiral position swallows a missing-file
    // error as "nothing catalogued here", so without this a wrong -d/-D reads
    // nothing everywhere and surfaces as InsufficientQuads - exit 1, "no
    // solution" - when the image is fine and the database is the problem.
    if !crate::catalog::catalog_present(&params.db_path, &params.db_name) {
        return Err(ArcsecError::CatalogNotFound(params.db_path.clone()));
    }

    // --- Phase A: star detection ---
    let bg = get_background(img, params.max_stars);
    log::info!("Start finding stars");
    let (stars, stars_raw) = find_stars_with_background(
        img,
        &bg,
        params.hfd_min,
        params.max_stars,
        img.width,
        img.height,
    );
    log::info!(
        "{} stars found of the requested {}. Background value is {:.0}. \
         Detection level used {:.0} above background. Star level is {:.0} above background. \
         Noise level is {:.0}",
        stars_raw,
        params.max_stars,
        bg.mean,
        bg.star_level,
        bg.star_level,
        bg.noise,
    );
    if stars_raw > params.max_stars {
        log::info!("Selecting the {} brightest stars only.", params.max_stars);
    }

    // Detection is not trimmed. Stars beyond `-s` are faint enough to be absent
    // from the catalog, which once corrupted 3-NN quads badly enough to justify
    // dropping all but the brightest half; quad redundancy and star-level
    // verification absorb that now, and halving the list halved the quad count.
    // The catalog still reads the full requested depth.

    let nrstars_image = stars.len();
    if nrstars_image < 5 {
        return Err(ArcsecError::InsufficientStars {
            found: nrstars_image,
            required: 5,
        });
    }

    // --- Phase B: image pattern building ---
    let img_quads = build_quads(&stars, nrstars_image);
    let nr_quads = img_quads.len();

    let img_tris = if params.method == SolveMethod::Tetra {
        build_triangles(&stars)
    } else {
        crate::quads::TriangleList::default()
    };

    let patterns_empty = match params.method {
        SolveMethod::Quads => nr_quads == 0,
        SolveMethod::Tetra => img_tris.is_empty(),
    };
    if patterns_empty {
        return Err(ArcsecError::InsufficientQuads {
            found: 0,
            required: 3,
        });
    }

    let min_quads: usize = 3 + nrstars_image / 140;

    let oversize: f64 = match params.speed {
        SearchSpeed::Auto if nrstars_image < 35 => 2.0,
        SearchSpeed::Auto if nrstars_image > 140 => 1.0,
        SearchSpeed::Auto => 2.0 * (35.0 / nrstars_image as f64).sqrt(),
        // As ASTAP, never more than one database tile: a larger window could reach
        // past the neighbouring tile, which the tile lookup does not cover.
        SearchSpeed::Slow => {
            let max_fov_deg = match crate::catalog::detect_layout(&params.db_path, &params.db_name)
            {
                CatalogLayout::Areas1476 => 5.142_857_143_f64,
                CatalogLayout::Areas290 => 9.53,
                CatalogLayout::AllSky001 => 180.0,
            };
            2.0_f64.min(max_fov_deg.to_radians() / params.fov).max(1.0)
        }
    };

    // Use the full catalog depth regardless of how many image stars we trimmed.
    let nrstars_required = (params.max_stars as f64 * oversize * oversize).round() as usize;
    let step_size = params.fov;
    let fov_deg = step_size.to_degrees();
    let max_distance = (params.search_radius / step_size + 2.0) as i32;

    log::info!(
        "{} stars, {} quads selected in the image. {} database stars, {} database quads required \
         for the {:.2}d square search window. Step size {:.2}d. Oversize {:.2}",
        nrstars_image,
        nr_quads,
        nrstars_required,
        nrstars_required,
        fov_deg * oversize,
        fov_deg,
        oversize,
    );

    // --- Phase C: spiral search ---
    //
    // Spiral positions are independent, so they are evaluated a batch at a time
    // across a thread pool. Semantics are unchanged from the serial search: within a
    // batch the lowest spiral index wins, and batches are processed in order, so the
    // position returned is exactly the one the serial loop would have returned. The
    // only cost is evaluating the rest of a batch after its first success.
    let ctx = SpiralCtx {
        params,
        img,
        stars: &stars,
        img_quads: &img_quads,
        img_tris: &img_tris,
        nrstars_image,
        nrstars_required,
        oversize,
        min_quads,
        step_size,
    };

    let n_threads = if params.threads > 0 {
        params.threads
    } else {
        crate::max_threads()
    }
    .clamp(1, 64);

    let positions: Vec<(i32, i32)> = SpiralSearch::new(max_distance).collect();
    let mut step_distances: Vec<f64> = Vec::new();

    let mut winner: Option<PositionOutcome> = None;
    let mut start_idx = 0usize;
    while start_idx < positions.len() && winner.is_none() {
        // The first position is the hint itself and usually solves outright, so try it
        // on its own: spawning a pool for it would cost more than it saves.
        let batch_len = if start_idx == 0 {
            1
        } else {
            n_threads.min(positions.len() - start_idx)
        };
        let batch = &positions[start_idx..start_idx + batch_len];

        let tries: Vec<PositionTry> = if n_threads == 1 || batch.len() == 1 {
            batch
                .iter()
                .enumerate()
                .map(|(k, &(sx, sy))| try_position(&ctx, start_idx + k, sx, sy))
                .collect()
        } else {
            std::thread::scope(|scope| {
                let handles: Vec<_> = batch
                    .iter()
                    .enumerate()
                    .map(|(k, &(sx, sy))| {
                        let ctx = &ctx;
                        scope.spawn(move || try_position(ctx, start_idx + k, sx, sy))
                    })
                    .collect();
                handles
                    .into_iter()
                    // A dead worker must not read as "nothing matched here":
                    // the spiral would move on and the solve would fail for a
                    // reason with no trace anywhere.
                    .map(|h| h.join().unwrap_or_else(|e| std::panic::resume_unwind(e)))
                    .collect()
            })
        };

        for t in tries {
            if let Some(d) = t.sep_deg {
                step_distances.push(d);
            }
            if let Some(o) = t.outcome
                && winner.as_ref().is_none_or(|w| o.idx < w.idx)
            {
                winner = Some(o);
            }
        }

        start_idx += batch_len;
    }

    if let Some(o) = winner {
        log::info!(
            "{} of {} patterns selected matching within {:.3} tolerance.",
            o.n_matched,
            o.n_raw,
            params.quad_tolerance,
        );

        let v = o.verified;
        let mut wcs = derive_wcs(o.ra_db, o.dec_db, &v.plate, img.width, img.height);
        // The verified pairs, on the original image's pixel grid: a binned pixel
        // centre at 0-based `x` is at `(x + 0.5) * b + 0.5` in unbinned FITS pixels.
        let b = params.binning.max(1) as f64;
        wcs.matched_stars = v
            .img_pos
            .iter()
            .zip(&v.cat_pos)
            .map(|(&(x, y), &(sx, sy))| {
                let (ra, dec) = standard_equatorial(o.ra_db, o.dec_db, sx, sy, 1.0);
                MatchedStar {
                    x: (x + 0.5) * b + 0.5,
                    y: (y + 0.5) * b + 0.5,
                    ra,
                    dec,
                }
            })
            .collect();
        if params.binning > 1 {
            let b = params.binning as f64;
            wcs.crpix1 = (wcs.crpix1 - 0.5) * b + 0.5;
            wcs.crpix2 = (wcs.crpix2 - 0.5) * b + 0.5;
            wcs.cd1_1 /= b;
            wcs.cd1_2 /= b;
            wcs.cd2_1 /= b;
            wcs.cd2_2 /= b;
            wcs.cdelt1 /= b;
            wcs.cdelt2 /= b;
        }
        wcs.residual_rms = v.rms;
        wcs.stars_matched = v.n();
        wcs.raw_matches = o.n_raw;
        wcs.plate = v.plate;
        wcs.mag_limit = o.mag_limit;
        wcs.search_dist_deg = o.sep_deg;
        wcs.step_distances = step_distances;
        return Ok(wcs);
    }

    Err(ArcsecError::InsufficientQuads {
        found: 0,
        required: min_quads,
    })
}

/// Format RA (radians) as `astap_cli` prints it: `"HH: MM  SS.S"`, each field at
/// least two digits (ASTAP's `prepare_ra(ra, ': ')`).
#[must_use]
pub fn format_ra(ra_rad: f64) -> String {
    // Round once, at the printed precision, and only then split into fields.
    // Splitting first and letting `{:.1}` round the seconds printed 59.96 s as
    // "60.0" without carrying into the minutes (and 23:59:59.96 as "23: 59  60.0").
    // ASTAP carries, but prints 23:59:59.96 as "24: 00  00.0"; this wraps to 00h.
    const TENTHS_PER_DAY: f64 = 24.0 * 36_000.0;
    let ra_tenths = ((ra_rad.to_degrees() / 15.0 * 36_000.0)
        .round()
        .rem_euclid(TENTHS_PER_DAY)) as u64;
    let h = ra_tenths / 36_000;
    let m = ra_tenths / 600 % 60;
    let s = ra_tenths % 600 / 10;
    let tenths = ra_tenths % 10;
    format!("{h:02}: {m:02}  {s:02}.{tenths}")
}

/// Format Dec (radians) as `astap_cli` prints it: `"±DDd MM  SS"`, each field at
/// least two digits (ASTAP's `prepare_dec(dec, 'd ')`).
#[must_use]
pub fn format_dec(dec_rad: f64) -> String {
    let dec_deg = dec_rad.to_degrees();
    let sign = if dec_deg < 0.0 { '-' } else { '+' };
    let dec_secs = (dec_deg.abs() * 3600.0).round() as u64;
    let dd = dec_secs / 3600;
    let dm = dec_secs / 60 % 60;
    let ds = dec_secs % 60;
    format!("{sign}{dd:02}d {dm:02}  {ds:02}")
}

/// Format RA and Dec (radians) as `astap_cli`'s `Solution found:` line does:
/// `"HH: MM  SS.S ±DDd MM  SS"`. (Its `Start position:` line puts a comma between
/// the two; see [`format_ra`] and [`format_dec`].)
#[must_use]
pub fn format_radec(ra_rad: f64, dec_rad: f64) -> String {
    format!("{} {}", format_ra(ra_rad), format_dec(dec_rad))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::math::coords::{ang_sep, standard_equatorial};
    use crate::test_support::{
        Rng, SkySpec, TempDir, TruthWcs, random_sky, render, write_001_db, write_290_db,
        write_1476_db,
    };
    use crate::types::{ImageBuffer, PlateConstants};
    use crate::wcs::output::derive_wcs;
    use core::f64::consts::PI;

    fn deg(d: f64) -> f64 {
        d * PI / 180.0
    }

    fn make_test_scene(
        n_stars: usize,
        ra_center: f64,
        dec_center: f64,
        cdelt_arcsec: f64,
        width: usize,
        height: usize,
    ) -> (ImageBuffer, Vec<(f64, f64)>, PlateConstants) {
        let mut data = vec![100.0f32; width * height];
        let mut catalog_sky: Vec<(f64, f64)> = Vec::new();
        let stars_per_row = (n_stars as f64).sqrt().ceil() as usize;
        let spacing = 40.0;
        let cx = (width as f64 - 1.0) / 2.0;
        let cy = (height as f64 - 1.0) / 2.0;
        let a = cdelt_arcsec;
        let c = -a * cx;
        let e = cdelt_arcsec;
        let f_offset = -e * cy;
        let plate = PlateConstants {
            a,
            b: 0.0,
            c,
            d: 0.0,
            e,
            f: f_offset,
        };
        let mut count = 0;
        'outer: for row in 0..stars_per_row {
            for col in 0..stars_per_row {
                if count >= n_stars {
                    break 'outer;
                }
                let px = 20.0 + col as f64 * spacing;
                let py = 20.0 + row as f64 * spacing;
                if px >= width as f64 - 20.0 || py >= height as f64 - 20.0 {
                    continue;
                }
                let x_std = a * px + c;
                let y_std = e * py + f_offset;
                let (ra, dec) = standard_equatorial(ra_center, dec_center, x_std, y_std, 1.0);
                catalog_sky.push((ra, dec));
                let sigma = 2.0;
                let amp = 30000.0f32;
                for dy in -8i32..=8 {
                    for dx in -8i32..=8 {
                        let x = (px as i32 + dx) as usize;
                        let y = (py as i32 + dy) as usize;
                        if x < width && y < height {
                            let r2 = (dx * dx + dy * dy) as f64 / (2.0 * sigma * sigma);
                            data[y * width + x] += amp * (-r2).exp() as f32;
                        }
                    }
                }
                count += 1;
            }
        }
        let img = ImageBuffer {
            data,
            width,
            height,
        };
        (img, catalog_sky, plate)
    }

    #[test]
    fn derive_wcs_recovers_position() {
        let ra_center = deg(45.0);
        let dec_center = deg(30.0);
        let (img, _cat, plate) = make_test_scene(16, ra_center, dec_center, 2.0, 300, 300);
        let wcs = derive_wcs(ra_center, dec_center, &plate, img.width, img.height);
        let sep_arcsec = ang_sep(wcs.ra0, wcs.dec0, ra_center, dec_center) * (180.0 / PI * 3600.0);
        assert!(sep_arcsec < 0.5, "centre offset = {sep_arcsec} arcsec");
    }

    #[test]
    fn spiral_covers_origin_first() {
        assert_eq!(SpiralSearch::new(5).next(), Some((0, 0)));
    }

    #[test]
    fn oversize_formula_limits() {
        for n in [10, 35, 70, 140, 200] {
            let ov: f64 = if n < 35 {
                2.0
            } else if n > 140 {
                1.0
            } else {
                2.0 * (35.0 / n as f64).sqrt()
            };
            assert!((1.0..=2.0).contains(&ov), "oversize={ov} for n={n}");
        }
    }

    #[test]
    fn format_radec_carries_rounded_seconds() {
        // 1h 59m 59.97s must round up to 2h 00m 00.0s, not print "60.0" seconds.
        let ra = deg((1.0 + 59.0 / 60.0 + 59.97 / 3600.0) * 15.0);
        // +10° 59' 59.7" rounds to +11° 00' 00".
        let dec = deg(10.0 + 59.0 / 60.0 + 59.7 / 3600.0);
        assert_eq!(format_radec(ra, dec), "02: 00  00.0 +11d 00  00");
        // RA just short of 24h wraps to 0h.
        let s = format_radec(deg(359.999_999_9), deg(-0.5));
        assert_eq!(s, "00: 00  00.0 -00d 30  00");
        // An ordinary value is unchanged by the rewrite.
        assert_eq!(
            format_radec(deg((5.0 + 35.0 / 60.0 + 17.3 / 3600.0) * 15.0), deg(-5.39)),
            "05: 35  17.3 -05d 23  24"
        );
    }

    /// Byte-for-byte what `astap_cli` prints (checked against 2026.07.30):
    /// `Start position: 04: 20  00.0, +35d 00  00` and
    /// `Solution found: 04: 20  00.0 +35d 00  00`. Every field is at least two
    /// digits wide, as ASTAP's `LeadingZero` makes it.
    #[test]
    fn ra_and_dec_are_formatted_as_astap_cli_prints_them() {
        let ra = deg(65.0); // 4h 20m
        let dec = deg(35.0);
        assert_eq!(format_ra(ra), "04: 20  00.0");
        assert_eq!(format_dec(dec), "+35d 00  00");
        assert_eq!(format_radec(ra, dec), "04: 20  00.0 +35d 00  00");
        assert_eq!(
            format_radec(
                deg((13.0 + 7.0 / 60.0 + 9.25 / 3600.0) * 15.0),
                -deg(89.0 + 1.0 / 60.0 + 2.0 / 3600.0)
            ),
            "13: 07  09.3 -89d 01  02"
        );
        assert_eq!(format_dec(deg(-0.0001)), "-00d 00  00");
    }

    #[test]
    fn solve_image_rejects_a_non_positive_fov() {
        let img = ImageBuffer::new(64, 64);
        let params = SolveParams {
            ra_hint: 0.0,
            dec_hint: 0.0,
            fov: 0.0,
            search_radius: 0.1,
            quad_tolerance: 0.007,
            hfd_min: 1.5,
            max_stars: 500,
            db_path: std::path::PathBuf::from("/nonexistent"),
            db_name: "d50".into(),
            binning: 1,
            method: SolveMethod::Quads,
            threads: 1,
            speed: SearchSpeed::Auto,
        };
        assert!(matches!(
            solve_image(&img, &params),
            Err(ArcsecError::InvalidParameter(_))
        ));
    }

    // ── Plate-fit helpers ─────────────────────────────────────────────────────

    /// A known similarity transform (pixels → catalogue arcsec), with a flip.
    fn known_plate() -> PlateConstants {
        let (s, r) = (3.2_f64, 0.61_f64);
        PlateConstants {
            a: -s * r.cos(),
            b: s * r.sin(),
            c: 640.0,
            d: s * r.sin(),
            e: s * r.cos(),
            f: -512.0,
        }
    }

    fn apply(p: &PlateConstants, (x, y): (f64, f64)) -> (f64, f64) {
        (p.a * x + p.b * y + p.c, p.d * x + p.e * y + p.f)
    }

    fn plate_close(p: &PlateConstants, q: &PlateConstants, tol: f64) -> bool {
        [
            (p.a, q.a),
            (p.b, q.b),
            (p.c, q.c),
            (p.d, q.d),
            (p.e, q.e),
            (p.f, q.f),
        ]
        .iter()
        .all(|(u, v)| (u - v).abs() <= tol)
    }

    fn star_at(x: f64, y: f64) -> Star {
        Star {
            x,
            y,
            snr: 50.0,
            hfd: 2.5,
        }
    }

    /// 40 exact pairs under `known_plate`, then five pairs whose catalogue side is
    /// displaced by `outlier(k)`.
    fn pairs_with_outliers(outlier: impl Fn(usize, (f64, f64)) -> (f64, f64)) -> PairedPositions {
        let plate = known_plate();
        let mut rng = Rng::new(7);
        let mut img = Vec::new();
        let mut cat = Vec::new();
        for _ in 0..40 {
            let p = (rng.range(0.0, 500.0), rng.range(0.0, 500.0));
            img.push(p);
            cat.push(apply(&plate, p));
        }
        for k in 0..5 {
            let p = (rng.range(0.0, 500.0), rng.range(0.0, 500.0));
            img.push(p);
            cat.push(outlier(k, apply(&plate, p)));
        }
        (img, cat)
    }

    #[test]
    fn sigma_clip_pairs_rejects_outliers_and_keeps_the_rest() {
        // Five wrong pairings, each ~100" (30 px) from where the plate puts them.
        let (img, cat) = pairs_with_outliers(|k, (x, y)| {
            let a = k as f64 * 1.3;
            (x + 100.0 * a.cos(), y + 100.0 * a.sin())
        });
        let (ci, cc) = sigma_clip_pairs(img, cat, 3.0, 3);
        assert_eq!(ci.len(), 40, "all and only the true pairs survive");
        let fit = solve_plate_constants(&ci, &cc).unwrap();
        assert!(plate_close(&fit, &known_plate(), 1e-6), "{fit:?}");
    }

    /// `sigma_clip_pairs` gives up as soon as a fit fails, and the first fit is made
    /// on the contaminated set. Five gross outliers in 45 pairs are enough to skew
    /// that fit past the 10% x/y scale check in `solve_plate_constants`
    /// (`BadSolution`, ratio 1.135 here), so nothing is clipped and all 45 come
    /// back. In `try_position` the Tetra path then refits the same contaminated set,
    /// fails the same check, and abandons a position whose 40 good pairs would
    /// have solved it. The clipper does not work in exactly the case it exists for.
    #[test]
    fn sigma_clip_pairs_rejects_gross_outliers() {
        let (img, cat) =
            pairs_with_outliers(|k, _| (1000.0 + 150.0 * k as f64, -900.0 + 70.0 * k as f64));
        assert!(matches!(
            solve_plate_constants(&img, &cat),
            Err(ArcsecError::BadSolution { .. })
        ));
        let (ci, _) = sigma_clip_pairs(img, cat, 3.0, 3);
        assert_eq!(ci.len(), 40, "the five gross outliers should be clipped");
    }

    #[test]
    fn sigma_clip_pairs_leaves_too_few_pairs_alone() {
        let img = vec![(0.0, 0.0), (1.0, 0.0)];
        let cat = vec![(5.0, 5.0), (9.0, 9.0)];
        let (ci, cc) = sigma_clip_pairs(img.clone(), cat.clone(), 3.0, 3);
        assert_eq!((ci, cc), (img, cat));
    }

    #[test]
    fn verify_and_refit_recovers_the_plate_from_a_rough_guess() {
        let truth = known_plate();
        let mut rng = Rng::new(11);
        let mut img_stars = Vec::new();
        let mut cat_stars = Vec::new();
        for _ in 0..60 {
            let (x, y) = (rng.range(5.0, 395.0), rng.range(5.0, 295.0));
            img_stars.push(star_at(x, y));
            let (cx, cy) = apply(&truth, (x, y));
            cat_stars.push(star_at(cx, cy));
        }
        // Catalogue stars that fall outside the frame must be ignored, not paired.
        for k in 0..20 {
            let (cx, cy) = apply(&truth, (-300.0 - 10.0 * k as f64, 900.0));
            cat_stars.push(star_at(cx, cy));
        }
        // Start 2 px and a little rotation away from the truth.
        let mut rough = truth.clone();
        rough.c += 2.0 * truth.a;
        rough.f += 2.0 * truth.e;
        rough.b += 0.01;
        let v = verify_and_refit(&StarList(img_stars), &StarList(cat_stars), &rough, 400, 300)
            .expect("a correct plate must verify");
        assert_eq!(v.n(), 60);
        assert_eq!(v.cat_pos.len(), 60);
        assert!(v.rms < 1e-6, "rms {}", v.rms);
        assert!(plate_close(&v.plate, &truth, 1e-6), "{:?}", v.plate);
        // Each pair is a star and its own catalogue entry.
        for (&(x, y), &(cx, cy)) in v.img_pos.iter().zip(&v.cat_pos) {
            let (px, py) = apply(&truth, (x, y));
            assert!((px - cx).hypot(py - cy) < 1e-6);
        }
    }

    #[test]
    fn verify_and_refit_rejects_too_few_or_clustered_matches() {
        let truth = known_plate();
        let mut rng = Rng::new(12);
        let build = |pts: &[(f64, f64)]| {
            let img = StarList(pts.iter().map(|&(x, y)| star_at(x, y)).collect());
            let cat = StarList(
                pts.iter()
                    .map(|&p| apply(&truth, p))
                    .map(|(x, y)| star_at(x, y))
                    .collect(),
            );
            (img, cat)
        };

        // 20 well-spread stars: fewer than MIN_VERIFIED_STARS.
        let few: Vec<_> = (0..20)
            .map(|_| (rng.range(0.0, 400.0), rng.range(0.0, 300.0)))
            .collect();
        let (img, cat) = build(&few);
        assert!(verify_and_refit(&img, &cat, &truth, 400, 300).is_none());

        // 80 stars, all in one 40-pixel corner: rotation is unconstrained.
        let clustered: Vec<_> = (0..80)
            .map(|_| (rng.range(0.0, 40.0), rng.range(0.0, 40.0)))
            .collect();
        let (img, cat) = build(&clustered);
        assert!(verify_and_refit(&img, &cat, &truth, 400, 300).is_none());

        // The same 80 spread over the frame pass.
        let spread: Vec<_> = (0..80)
            .map(|_| (rng.range(0.0, 400.0), rng.range(0.0, 300.0)))
            .collect();
        let (img, cat) = build(&spread);
        assert!(verify_and_refit(&img, &cat, &truth, 400, 300).is_some());

        // Degenerate inputs.
        let empty = StarList::default();
        assert!(verify_and_refit(&empty, &cat, &truth, 400, 300).is_none());
        let mut singular = truth.clone();
        singular.a = 0.0;
        singular.b = 0.0;
        assert!(verify_and_refit(&img, &cat, &singular, 400, 300).is_none());
    }

    // ── End-to-end solves against synthetic catalogues ────────────────────────

    #[derive(Clone, Copy)]
    enum Db {
        Areas1476,
        Areas290,
        AllSky001,
    }

    /// A rendered field and the database it was drawn from.
    struct Scene {
        dir: TempDir,
        img: ImageBuffer,
        truth: TruthWcs,
    }

    /// Render ~`n_in_frame` stars through `truth` and write the surrounding sky
    /// (six fields wide, so offset hints still find their stars) as a database.
    fn scene(truth: TruthWcs, db: Db, n_in_frame: usize, seed: u64) -> Scene {
        let mut rng = Rng::new(seed);
        let scale_deg = truth.cd[1].hypot(truth.cd[3]);
        let (w_deg, h_deg) = (
            truth.width as f64 * scale_deg,
            truth.height as f64 * scale_deg,
        );
        let side = 6.0 * w_deg.max(h_deg);
        let sky = random_sky(
            &mut rng,
            &SkySpec {
                ra0: truth.ra0,
                dec0: truth.dec0,
                side_deg: side,
                n: (n_in_frame as f64 * side * side / (w_deg * h_deg)) as usize,
                min_sep_deg: 12.0 * scale_deg,
                mag_lo: 10.0,
                mag_hi: 14.5,
            },
        );
        // A PSF of ~1.3 px on a 5"/px frame, scaled so binned frames stay sampled.
        let sigma = 1.3 * 5.0 / (scale_deg * 3600.0);
        let img = render(
            &truth,
            &sky,
            sigma.max(1.3),
            1000.0,
            8.0,
            30_000.0,
            &mut rng,
        );
        let dir = TempDir::new("solve");
        match db {
            Db::Areas1476 => write_1476_db(dir.path(), "t50", &sky),
            Db::Areas290 => write_290_db(dir.path(), "t50", &sky),
            Db::AllSky001 => write_001_db(dir.path(), "t50", &sky),
        }
        Scene { dir, img, truth }
    }

    fn params_for(s: &Scene, ra_hint: f64, dec_hint: f64) -> SolveParams {
        SolveParams {
            ra_hint,
            dec_hint,
            fov: (s.truth.height as f64 * s.truth.cd[1].hypot(s.truth.cd[3])).to_radians(),
            search_radius: deg(2.0),
            quad_tolerance: 0.007,
            hfd_min: 1.5,
            max_stars: 500,
            db_path: s.dir.path().to_path_buf(),
            db_name: "t50".into(),
            binning: 1,
            method: SolveMethod::Quads,
            threads: 1,
            speed: SearchSpeed::Auto,
        }
    }

    fn assert_solved(s: &Scene, wcs: &WcsSolution, tol_arcsec: f64) {
        let err = s.truth.max_error_arcsec(wcs);
        assert!(
            err < tol_arcsec,
            "worst centre/corner error {err:.3}\" (matched {}, rms {:.3})",
            wcs.stars_matched,
            wcs.residual_rms
        );
        assert!(wcs.stars_matched >= MIN_VERIFIED_STARS);
        // Star-level residual under a third of a pixel.
        let scale_arcsec = s.truth.cd[1].hypot(s.truth.cd[3]) * 3600.0;
        assert!(
            wcs.residual_rms < 0.3 * scale_arcsec,
            "rms {}",
            wcs.residual_rms
        );
        assert!(wcs.raw_matches > 0);
        assert_matches_agree(wcs, 0.3, 1.0);
        assert!(wcs.mag_limit > 10.0 && wcs.mag_limit <= 14.5);
        assert!(
            wcs.cdelt1 < 0.0 && wcs.cdelt2 > 0.0,
            "CDELT sign convention"
        );
    }

    /// The verified pairs agree with the solution: RMS under `rms_px` pixels, and
    /// none further off than the final verification radius (`binning` pixels each).
    fn assert_matches_agree(wcs: &WcsSolution, rms_px: f64, binning: f64) {
        assert_eq!(wcs.matched_stars.len(), wcs.stars_matched);
        assert!(wcs.sip.is_none(), "solve_image never fits SIP");
        let tan = crate::wcs::TanWcs::from(wcs);
        let mut sq = 0.0;
        for m in &wcs.matched_stars {
            let (x, y) = tan.sky_to_pixel(m.ra, m.dec).unwrap();
            let d = (x - m.x).hypot(y - m.y);
            assert!(
                d < VERIFY_RADII[VERIFY_RADII.len() - 1] * binning,
                "pair at ({:.2},{:.2}) projects to ({x:.2},{y:.2})",
                m.x,
                m.y
            );
            sq += d * d;
        }
        let rms = (sq / wcs.matched_stars.len() as f64).sqrt();
        assert!(rms < rms_px, "pair rms {rms} px");
    }

    #[test]
    fn solves_a_1476_database_from_an_offset_hint() {
        let truth = TruthWcs::new(deg(84.3), deg(-5.2), 5.0, 23.0, false, 400, 320);
        let s = scene(truth, Db::Areas1476, 130, 1);
        // Hint roughly one field away in each axis: the spiral has to move.
        let mut p = params_for(&s, deg(84.3 + 0.6), deg(-5.2 - 0.45));
        p.threads = 4;
        let wcs = solve_image(&s.img, &p).expect("solve");
        assert_solved(&s, &wcs, 1.0);
        assert!(wcs.search_dist_deg > 0.1, "solved at the hint itself?");
        assert!(wcs.step_distances.len() > 1);
        // The pixel scale and rotation come back too.
        assert!((wcs.cdelt2 * 3600.0 - 5.0).abs() < 0.01, "{}", wcs.cdelt2);
        assert!((wcs.crota2 - 23.0).abs() < 0.05, "crota2 {}", wcs.crota2);
    }

    #[test]
    fn solves_a_mirrored_image_on_a_290_database() {
        let truth = TruthWcs::new(deg(201.0), deg(47.5), 6.0, 160.0, true, 360, 360);
        let s = scene(truth, Db::Areas290, 120, 2);
        let wcs = solve_image(&s.img, &params_for(&s, truth.ra0, truth.dec0)).expect("solve");
        assert_solved(&s, &wcs, 1.0);
        assert!(wcs.search_dist_deg < 1e-9, "should solve at the hint");
        // A mirrored image has det(CD) > 0.
        assert!(wcs.cd1_1 * wcs.cd2_2 - wcs.cd1_2 * wcs.cd2_1 > 0.0);
    }

    #[test]
    fn solves_across_ra_zero_with_an_all_sky_001_database() {
        // The field straddles RA 0h, so its catalogue stars sit either side of 2π.
        let truth = TruthWcs::new(deg(0.05), deg(21.0), 5.0, -70.0, false, 360, 300);
        let s = scene(truth, Db::AllSky001, 120, 3);
        let wcs = solve_image(&s.img, &params_for(&s, truth.ra0, truth.dec0)).expect("solve");
        assert_solved(&s, &wcs, 1.0);
    }

    #[test]
    fn solves_across_ra_zero_with_a_1476_database() {
        let truth = TruthWcs::new(deg(359.97), deg(-33.0), 5.0, 95.0, false, 360, 300);
        let s = scene(truth, Db::Areas1476, 120, 4);
        let wcs = solve_image(&s.img, &params_for(&s, truth.ra0, truth.dec0)).expect("solve");
        assert_solved(&s, &wcs, 1.0);
    }

    #[test]
    fn solves_a_field_near_the_celestial_pole() {
        let truth = TruthWcs::new(deg(40.0), deg(88.9), 5.0, 10.0, false, 360, 300);
        let s = scene(truth, Db::Areas1476, 120, 5);
        let wcs = solve_image(&s.img, &params_for(&s, truth.ra0, truth.dec0)).expect("solve");
        assert_solved(&s, &wcs, 1.0);
    }

    /// The plate constants are fitted in the tangent plane of the spiral position
    /// that matched (`ra_db`, `dec_db`), but `derive_wcs` then moves CRVAL to the
    /// image centre and keeps the CD matrix unchanged, as though the two tangent
    /// planes were the same. They are not, and the error grows linearly with the
    /// distance between the matched spiral position and the true field centre.
    ///
    /// Measured on this 1.5° × 1.25° field at 15"/px: worst-corner error 0.35" with
    /// the hint on the centre, 7.5" at 0.2° off, 14.5" at 0.4°, 21" at 0.6° (1.4 px),
    /// while the star-level RMS stays at 0.7-0.85" throughout — the verification
    /// cannot see it, because it runs in the same (offset) tangent plane. Spiral
    /// positions land up to half a step (half a field) from the truth, and a blind
    /// estimate can be a whole field off, so this is well inside normal use. 5" is
    /// the corner error `scripts/benchmark.py` counts as a false positive.
    #[test]
    fn accuracy_does_not_depend_on_the_hint_offset() {
        let truth = TruthWcs::new(deg(150.0), deg(30.0), 15.0, 20.0, false, 360, 300);
        let s = scene(truth, Db::Areas1476, 120, 21);
        let off = 0.4;
        let p = params_for(&s, deg(150.0 + off / deg(30.0).cos()), deg(30.0 + off));
        let wcs = solve_image(&s.img, &p).expect("solve");
        assert!(wcs.search_dist_deg < 1e-9, "solved at the hint");
        let err = s.truth.max_error_arcsec(&wcs);
        assert!(
            err < 5.0,
            "worst corner error {err:.2}\" with a {off}° hint offset"
        );
    }

    #[test]
    fn solves_with_the_tetra_method() {
        let truth = TruthWcs::new(deg(150.0), deg(2.0), 5.0, 45.0, false, 360, 300);
        let s = scene(truth, Db::Areas1476, 110, 6);
        let mut p = params_for(&s, truth.ra0, truth.dec0);
        p.method = SolveMethod::Tetra;
        let wcs = solve_image(&s.img, &p).expect("solve");
        assert_solved(&s, &wcs, 1.0);
    }

    #[test]
    fn slow_speed_solves_from_an_offset_hint() {
        let truth = TruthWcs::new(deg(84.3), deg(-5.2), 5.0, 23.0, false, 400, 320);
        let s = scene(truth, Db::Areas1476, 130, 1);
        let mut p = params_for(&s, deg(84.3 + 0.6), deg(-5.2 - 0.45));
        p.speed = SearchSpeed::Slow;
        let wcs = solve_image(&s.img, &p).expect("solve");
        assert_solved(&s, &wcs, 1.0);
    }

    #[test]
    fn binned_solve_is_reported_on_the_unbinned_pixel_grid() {
        // Render at full resolution, then solve the 2×2-binned frame.
        let truth = TruthWcs::new(deg(10.0), deg(40.0), 2.5, 30.0, false, 720, 600);
        let s = scene(truth, Db::Areas1476, 120, 7);
        let binned = s.img.bin_image(2);
        assert_eq!((binned.width, binned.height), (360, 300));
        let mut p = params_for(&s, truth.ra0, truth.dec0);
        p.binning = 2;
        let wcs = solve_image(&binned, &p).expect("solve");
        // crpix is the centre of the unbinned frame, and the scale is unbinned.
        assert!((wcs.crpix1 - 360.5).abs() < 1e-9, "crpix1 {}", wcs.crpix1);
        assert!((wcs.crpix2 - 300.5).abs() < 1e-9, "crpix2 {}", wcs.crpix2);
        assert!((wcs.cdelt2 * 3600.0 - 2.5).abs() < 0.01, "{}", wcs.cdelt2);
        let err = s.truth.max_error_arcsec(&wcs);
        assert!(err < 2.0, "worst corner error {err:.3}\"");
        // The pairs are on the unbinned grid too: one binned pixel is two of these.
        assert_matches_agree(&wcs, 0.6, 2.0);
    }

    #[test]
    fn a_field_absent_from_the_catalogue_does_not_solve() {
        // The image shows one random sky, the database holds a different one at the
        // same place: nothing may verify, however many quads happen to match.
        let truth = TruthWcs::new(deg(120.0), deg(-40.0), 5.0, 0.0, false, 360, 300);
        let s = scene(truth, Db::Areas1476, 120, 8);
        let decoy = TempDir::new("decoy");
        let mut rng = Rng::new(99);
        let other = random_sky(
            &mut rng,
            &SkySpec {
                ra0: truth.ra0,
                dec0: truth.dec0,
                side_deg: 3.0,
                n: 4000,
                min_sep_deg: 0.015,
                mag_lo: 10.0,
                mag_hi: 14.5,
            },
        );
        write_1476_db(decoy.path(), "t50", &other);
        let mut p = params_for(&s, truth.ra0, truth.dec0);
        p.db_path = decoy.path().to_path_buf();
        p.search_radius = deg(0.5);
        match solve_image(&s.img, &p) {
            Err(ArcsecError::InsufficientQuads { found: 0, required }) => {
                assert!(required >= 3);
            }
            other => panic!("expected InsufficientQuads, got {other:?}"),
        }
    }

    #[test]
    fn a_corrupt_catalogue_tile_is_skipped_not_fatal() {
        let truth = TruthWcs::new(deg(120.0), deg(-40.0), 5.0, 0.0, false, 360, 300);
        let s = scene(truth, Db::Areas1476, 120, 9);
        // Declare an unsupported record size in every tile.
        for entry in std::fs::read_dir(s.dir.path()).unwrap() {
            let path = entry.unwrap().path();
            let mut bytes = std::fs::read(&path).unwrap();
            bytes[109] = 7;
            std::fs::write(&path, bytes).unwrap();
        }
        let mut p = params_for(&s, truth.ra0, truth.dec0);
        p.search_radius = 0.0;
        assert!(matches!(
            solve_image(&s.img, &p),
            Err(ArcsecError::InsufficientQuads { .. })
        ));
    }

    #[test]
    fn a_blank_frame_reports_insufficient_stars() {
        let dir = TempDir::new("blank");
        write_1476_db(dir.path(), "t50", &[]);
        let mut rng = Rng::new(3);
        let img = ImageBuffer {
            data: (0..200 * 200)
                .map(|_| (1000.0 + 5.0 * rng.gauss()) as f32)
                .collect(),
            width: 200,
            height: 200,
        };
        let p = SolveParams {
            ra_hint: 0.0,
            dec_hint: 0.0,
            fov: deg(0.3),
            search_radius: deg(1.0),
            quad_tolerance: 0.007,
            hfd_min: 1.5,
            max_stars: 500,
            db_path: dir.path().to_path_buf(),
            db_name: "t50".into(),
            binning: 1,
            method: SolveMethod::Quads,
            threads: 1,
            speed: SearchSpeed::Auto,
        };
        match solve_image(&img, &p) {
            Err(ArcsecError::InsufficientStars { found, required: 5 }) => assert!(found < 5),
            other => panic!("expected InsufficientStars, got {other:?}"),
        }
    }

    #[test]
    fn a_missing_database_is_reported_before_any_detection() {
        let dir = TempDir::new("nodb");
        let p = SolveParams {
            ra_hint: 0.0,
            dec_hint: 0.0,
            fov: deg(1.0),
            search_radius: deg(1.0),
            quad_tolerance: 0.007,
            hfd_min: 1.5,
            max_stars: 500,
            db_path: dir.path().to_path_buf(),
            db_name: "d50".into(),
            binning: 1,
            method: SolveMethod::Quads,
            threads: 1,
            speed: SearchSpeed::Auto,
        };
        match solve_image(&ImageBuffer::new(64, 64), &p) {
            Err(ArcsecError::CatalogNotFound(path)) => assert_eq!(path, dir.path()),
            other => panic!("expected CatalogNotFound, got {other:?}"),
        }
    }

    #[test]
    fn solve_image_rejects_a_bad_search_radius_or_fov() {
        let base = SolveParams {
            ra_hint: 0.0,
            dec_hint: 0.0,
            fov: deg(1.0),
            search_radius: 0.1,
            quad_tolerance: 0.007,
            hfd_min: 1.5,
            max_stars: 500,
            db_path: std::path::PathBuf::from("/nonexistent"),
            db_name: "d50".into(),
            binning: 1,
            method: SolveMethod::Quads,
            threads: 1,
            speed: SearchSpeed::Auto,
        };
        let img = ImageBuffer::new(64, 64);
        for (fov, radius) in [
            (f64::NAN, 0.1),
            (-1.0, 0.1),
            (f64::INFINITY, 0.1),
            (0.01, -0.1),
            (0.01, f64::NAN),
            (0.01, f64::INFINITY),
        ] {
            let p = SolveParams {
                fov,
                search_radius: radius,
                ..base.clone()
            };
            assert!(
                matches!(solve_image(&img, &p), Err(ArcsecError::InvalidParameter(_))),
                "fov {fov}, radius {radius}"
            );
        }
    }

    #[test]
    fn format_radec_roundtrip() {
        let s = format_radec(deg(160.875), deg(-59.524));
        assert!(s.contains("10:"), "RA hours: {s}");
        assert!(s.contains('-'), "dec sign: {s}");
    }
}