stdpipe.astrometry_quad module

Quad-hash astrometry solver (Python-only) with practical upgrades:

1. Multi-scale quad pools with canonical descriptors: - fully canonical quad descriptors (astrometry.net-style A/B convention),

   so a hash hit directly provides the point correspondence for the fit

   • quads are partitioned into baseline-length rank bins (e.g. 6 bins)

   • consistent binning between detection and reference sets via shared edges

   • matching is done only within the same bin, reducing ambiguity

2. Two-stage hypothesis scoring with accuracy enhancements: - Stage 1: weighted scoring on a small subset with multi-probe hashing - Stage 2: mutual nearest-neighbor matching with weighted scoring - Iterative affine re-matching to grow the match set - Progressive sigma-clipping in final WCS refinement

Inputs

obj : astropy.table.Table with ‘x’,’y’ and either ‘flux’ or ‘mag’ cat : astropy.table.Table with ‘ra’,’dec’,’mag’ wcs_init : astropy.wcs.WCS (rough)

Outputs

refined_wcs : astropy.wcs.WCS match : astropy.table.Table diagnostics : dict

Deps: numpy, scipy, astropy

stdpipe.astrometry_quad.mag_orderings_consistent(mags_a, mags_b, tol: float = 0.5) → bool

Check that pairwise brightness orderings of two corresponding point sets agree within a tolerance.

Only pairs whose magnitude difference exceeds tol in both sets are required to agree in sign, so noise-flipped orderings of nearly equal-brightness stars are not rejected. Insensitive to a constant magnitude offset (instrumental vs catalogue zero point).

stdpipe.astrometry_quad.estimate_similarity_2d(A: ndarray, B: ndarray, allow_reflection: bool = True) → Tuple[ndarray, ndarray, float]

stdpipe.astrometry_quad.apply_similarity(xy: ndarray, R: ndarray, t: ndarray, s: float) → ndarray

stdpipe.astrometry_quad.estimate_affine_2d(A: ndarray, B: ndarray) → Tuple[ndarray, ndarray]

Fit affine transform B = A @ M.T + t via least squares (6 DOF).

More flexible than similarity (4 DOF) - handles shear and non-square pixels. Returns (M, t) where M is (2,2) linear part and t is (2,) translation.

stdpipe.astrometry_quad.apply_affine(xy: ndarray, M: ndarray, t: ndarray) → ndarray

Apply affine transform: result = xy @ M.T + t

Raises ValueError if the result contains non-finite values.

stdpipe.astrometry_quad.make_local_quads(points: ndarray, k: int = 8) → List[Tuple[int, int, int, int]]

stdpipe.astrometry_quad.quad_descriptor(P: ndarray) → Tuple[ndarray, float]

Descriptor invariant to translation/rotation/scale. Returns (desc[4], baseline_len).

stdpipe.astrometry_quad.quad_descriptor_batch(points: ndarray, quad_indices: ndarray) → Tuple[ndarray, ndarray, ndarray]

Vectorized canonical quad descriptor calculation for multiple quads.

The descriptor is a pure geometric invariant of the 4-point set: the baseline endpoint A is chosen such that cx + dx <= 1 (astrometry.net convention) and C, D are ordered lexicographically, so the same physical quad yields the same descriptor and point order regardless of the order in which the indices are listed.

Args:

points: (N, 2) array of point coordinates quad_indices: (M, 4) array of quad indices into points

Returns:

descs: (M, 4) array of descriptors lens: (M,) array of baseline lengths order: (M, 4) array of point indices in canonical (A, B, C, D)

order (-1 for degenerate quads)

stdpipe.astrometry_quad.quantize_desc(desc: ndarray, eps: float) → Tuple[int, int, int, int]

stdpipe.astrometry_quad.reflect_desc(desc: ndarray) → ndarray

Reflect quad descriptor across its baseline (flip y) with lexicographic re-ordering.

stdpipe.astrometry_quad.multiprobe_desc_keys(desc: ndarray, eps: float, threshold: float = 0.3) → List[Tuple[int, int, int, int]]

Multi-probe hashing: return the primary quantized key plus keys for nearby bins where the descriptor is close to a bin boundary.

Instead of searching all (2r+1)^4 = 81 neighboring bins, this probes only the bins that the descriptor might fall into due to quantization noise. Typically returns 1-8 keys instead of 81.

Args:

desc: 4D descriptor array eps: quantization step size threshold: fraction of eps within which a boundary is considered “near”

stdpipe.astrometry_quad.compute_bin_edges(lengths: ndarray, n_bins: int) → ndarray | None

Compute quantile-based bin edges from a set of baseline lengths.

Returns array of (n_bins+1) edge values, or None if binning not possible.

stdpipe.astrometry_quad.baseline_rank_bins(lengths: ndarray, n_bins: int, edges: ndarray | None = None) → ndarray

Assign each length to a rank bin based on quantiles. If edges are provided, use those instead of computing from the data. This allows consistent binning between detection and reference sets.

class stdpipe.astrometry_quad.QuadEntry(idxs: 'Tuple[int, int, int, int]', mags: 'Optional[Tuple[float, float, float, float]]' = None)

Bases: object

Attributes:

mags

idxs: Tuple[int, int, int, int]

mags: Tuple[float, float, float, float] | None = None

stdpipe.astrometry_quad.build_quad_hash_multiscale(points: ndarray, mags: ndarray | None, k: int, eps: float, n_scale_bins: int, allow_reflection: bool = True, bin_edges: ndarray | None = None, quad_array: ndarray | None = None, descs: ndarray | None = None, lens: ndarray | None = None, order: ndarray | None = None) → Dict[Tuple[int, Tuple[int, int, int, int]], List[QuadEntry]]

Hash key: (scale_bin_id, quantized_descriptor)

Entries store point indices in the canonical (A, B, C, D) descriptor order, so a hash hit directly provides the point correspondence for the transform fit.

If bin_edges is provided, use those for baseline binning instead of computing quantiles from the data. This enables consistent binning between detection and reference sets.

class stdpipe.astrometry_quad.AstrometryConfig(n_det: 'int' = 150, n_ref: 'int' = 600, neighbor_k: 'int' = 10, eps_desc: 'float' = 0.02, n_scale_bins: 'int' = 6, allow_reflection: 'bool' = True, use_mag_signature: 'bool' = True, mag_tolerance: 'float' = 0.5, use_wcs_prior: 'bool' = True, scale_tolerance: 'float' = 0.3, enforce_parity: 'bool' = True, stage1_n_det: 'int' = 60, stage1_radius_arcsec: 'float' = 20.0, stage2_radius_arcsec: 'float' = 8.0, top_k_hypotheses: 'int' = 80, max_quads_tested: 'int' = 8000, max_candidates_per_bucket: 'int' = 60, sip_degree: 'int' = 3, refine_clip_sigma: 'float' = 4.0, refine_clip_sigma_start: 'float' = 5.0, refine_min_match_fraction: 'float' = 0.5, refine_max_iter: 'int' = 5, refine_rematch_iters: 'int' = 2, refine_use_all: 'bool' = False, refine_n_det: 'Optional[int]' = None, refine_n_ref: 'Optional[int]' = None, refine_match_radius_arcsec: 'Optional[float]' = None, auto_match_resolution: 'bool' = False, adaptive_n_retry: 'int' = 2, adaptive_min_inliers: 'int' = 12)

Bases: object

Attributes:

refine_match_radius_arcsec

refine_n_det

refine_n_ref

n_det: int = 150

n_ref: int = 600

neighbor_k: int = 10

eps_desc: float = 0.02

n_scale_bins: int = 6

allow_reflection: bool = True

use_mag_signature: bool = True

mag_tolerance: float = 0.5

use_wcs_prior: bool = True

scale_tolerance: float = 0.3

enforce_parity: bool = True

stage1_n_det: int = 60

stage1_radius_arcsec: float = 20.0

stage2_radius_arcsec: float = 8.0

top_k_hypotheses: int = 80

max_quads_tested: int = 8000

max_candidates_per_bucket: int = 60

sip_degree: int = 3

refine_clip_sigma: float = 4.0

refine_clip_sigma_start: float = 5.0

refine_min_match_fraction: float = 0.5

refine_max_iter: int = 5

refine_rematch_iters: int = 2

refine_use_all: bool = False

refine_n_det: int | None = None

refine_n_ref: int | None = None

refine_match_radius_arcsec: float | None = None

auto_match_resolution: bool = False

adaptive_n_retry: int = 2

adaptive_min_inliers: int = 12

class stdpipe.astrometry_quad.QuadHashAstrometry(config: AstrometryConfig | None = None)

Bases: object

Methods

refine

refine(obj: Table, cat: Table, wcs_init: WCS, fit_projection: WCS | None = None, pv_deg: int = 5) → Tuple[WCS, Table, dict]

stdpipe.astrometry_quad.refine_wcs_quadhash(obj: Table, cat: Table, wcs: WCS = None, header=None, sr: float = 0.002777777777777778, order: int = 2, projection: str = None, pv_deg: int = 5, cat_col_ra: str = 'RAJ2000', cat_col_dec: str = 'DEJ2000', cat_col_mag: str = 'rmag', obj_col_mag: str = None, obj_col_flux: str = 'flux', sn: float = None, n_det: int = 150, n_ref: int = 600, max_quads_tested: int = 8000, allow_reflection: bool = True, use_wcs_prior: bool = True, scale_tolerance: float = 0.3, enforce_parity: bool = True, refine_use_all: bool = False, refine_n_det: int | None = None, refine_n_ref: int | None = None, refine_match_radius_arcsec: float | None = None, get_header: bool = False, update: bool = False, verbose: bool = False) → WCS

Refine WCS using quad-hash pattern matching algorithm.

Pure Python implementation of astrometric refinement using quad-hash based pattern matching with no external dependencies (only numpy, scipy, astropy). Typically achieves sub-arcsecond accuracy.

Parameters:

objastropy.table.Table

List of objects on the frame, must contain at least x, y, and either flux or mag columns.

catastropy.table.Table

Reference astrometric catalogue with RA, Dec, and magnitude columns.

wcsastropy.wcs.WCS, optional

Initial WCS solution (rough estimate).

headerastropy.io.fits.Header, optional

FITS header containing the initial astrometric solution (alternative to wcs).

srfloat, optional

Matching radius in degrees for stage 2 matching.

orderint, optional

SIP polynomial order for the distortion solution (0–5). For TAN projection, controls the SIP polynomial degree. For ZPN projection, controls additional SIP corrections on top of the PV radial model (0 means ZPN-only, 2 is recommended for wide-field images). Ignored for other projections.

projectionstr or None, optional

Target projection type for the output WCS. If None (default), uses the same projection as the input WCS. Supported values:

• 'TAN' — gnomonic (TAN) with SIP distortion polynomials

• 'ZPN' — zenithal polynomial with PV radial terms (+ optional SIP for non-radial corrections). Best for wide-field images (FoV > 5°) as ZPN naturally handles radial distortion that TAN-SIP struggles with.

• Any other FITS WCS projection code supported by astropy ('STG', 'ARC', 'ZEA', 'SIN', etc.) — linear WCS fit only (no distortion polynomials).

The initial WCS (wcs) is always used for the pattern matching phase; the projection conversion only affects the final WCS fit.

pv_degint, optional

ZPN PV polynomial degree (PV2_0 … PV2_pv_deg). Only used when projection='ZPN'. Default is 5, which is sufficient for most optical systems.

cat_col_rastr, optional

Catalogue column name for Right Ascension. Default is 'RAJ2000'.

cat_col_decstr, optional

Catalogue column name for Declination. Default is 'DEJ2000'.

cat_col_magstr, optional

Catalogue column name for magnitude. Default is 'rmag'.

obj_col_magstr, optional

Object list column name for magnitude. Auto-detected if not provided.

obj_col_fluxstr, optional

Object list column name for flux. Used if magnitude column is absent. Default is 'flux'.

snfloat, optional

If provided, only objects with S/N exceeding this value will be used.

n_detint, optional

Number of brightest detections to use for matching. Default is 150.

n_refint, optional

Number of brightest catalogue stars to use for matching. Default is 600.

max_quads_testedint, optional

Maximum number of detection quads to test. Default is 8000.

allow_reflectionbool, optional

Allow reflected quad matches. Default is True.

use_wcs_priorbool, optional

Filter transform hypotheses by scale (close to 1) and parity consistency in pixel space. Default is True.

scale_tolerancefloat, optional

Fractional tolerance for scale filtering. Default is 0.30 (±30%).

enforce_paritybool, optional

Enforce parity consistency with the initial WCS. Default is True.

refine_use_allbool, optional

If True, expand the final refinement to a larger pool of objects/catalogue entries.

refine_n_detint or None, optional

Max number of detections in the expanded refinement pool. None means all.

refine_n_refint or None, optional

Max number of catalogue stars in the expanded refinement pool. None means all.

refine_match_radius_arcsecfloat or None, optional

Matching radius in arcsec for the expanded refinement pool. Defaults to the stage 2 radius.

get_headerbool, optional

If True, return the FITS header object instead of WCS solution.

updatebool, optional

If True, update object sky coordinates in-place using the refined WCS.

verbosebool or callable, optional

Whether to show verbose messages. May be boolean or a print-like callable.

Returns:

astropy.wcs.WCS or astropy.io.fits.Header or None

Refined WCS solution, or FITS header if get_header=True, or None on failure.

Examples

>>> from stdpipe.astrometry_quad import refine_wcs_quadhash
>>> wcs_refined = refine_wcs_quadhash(
...     obj, cat, wcs_init,
...     sr=10/3600,  # 10 arcsec matching radius
...     order=2,     # quadratic SIP distortion
...     sn=5,        # use only S/N > 5 objects
...     verbose=True
... )

For wide-field images (FoV > 5 degrees), ZPN projection with SIP corrections typically gives lower systematic residuals at field edges:

>>> wcs_refined = refine_wcs_quadhash(
...     obj, cat, wcs_init,
...     sr=30/3600,
...     order=2,           # SIP order on top of ZPN radial model
...     projection='ZPN',  # zenithal polynomial projection
...     pv_deg=5,          # PV radial polynomial degree
...     verbose=True
... )