stdpipe.astrometry_quad module

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

1. Multi-scale quad pools: - 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_signature(mags4: ndarray) → Tuple[int, int, int, int]

stdpipe.astrometry_quad.tan_project_deg(ra_deg: ndarray, dec_deg: ndarray, ra0_deg: float, dec0_deg: float) → Tuple[ndarray, ndarray]

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]

Vectorized quad descriptor calculation for multiple quads.

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

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.neighbors_bins(key: Tuple[int, int, int, int], r: int = 1) → Iterable[Tuple[int, int, int, int]]

Generate neighbor bins using optimized pre-computed offsets.

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]', desc: 'np.ndarray', bin_id: 'int', mags: 'Optional[Tuple[float, float, float, float]]' = None)

Bases: object

Attributes:

mags

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

desc: ndarray

bin_id: 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) → Dict[Tuple[int, Tuple[int, int, int, int]], List[QuadEntry]]

Hash key: (scale_bin_id, quantized_descriptor)

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.

Phase 2 Optimization: Uses vectorized quad descriptor calculation.

class stdpipe.astrometry_quad.AstrometryConfig(n_det: 'int' = 180, n_ref: 'int' = 380, neighbor_k: 'int' = 8, eps_desc: 'float' = 0.015, bin_neighbor_radius: 'int' = 1, n_scale_bins: 'int' = 6, allow_reflection: 'bool' = True, use_mag_signature: 'bool' = True, use_wcs_prior: 'bool' = True, scale_tolerance: 'float' = 0.3, enforce_parity: 'bool' = True, stage1_n_det: 'int' = 50, stage1_radius_arcsec: 'float' = 20.0, stage2_radius_arcsec: 'float' = 8.0, top_k_hypotheses: 'int' = 60, max_quads_tested: 'int' = 5000, max_candidates_per_bucket: 'int' = 60, use_multiprobe: 'bool' = True, 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' = True, 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 = 180

n_ref: int = 380

neighbor_k: int = 8

eps_desc: float = 0.015

bin_neighbor_radius: int = 1

n_scale_bins: int = 6

allow_reflection: bool = True

use_mag_signature: bool = True

use_wcs_prior: bool = True

scale_tolerance: float = 0.3

enforce_parity: bool = True

stage1_n_det: int = 50

stage1_radius_arcsec: float = 20.0

stage2_radius_arcsec: float = 8.0

top_k_hypotheses: int = 60

max_quads_tested: int = 5000

max_candidates_per_bucket: int = 60

use_multiprobe: bool = True

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 = True

adaptive_n_retry: int = 2

adaptive_min_inliers: int = 12

class stdpipe.astrometry_quad.QuadHashAstrometry(config: AstrometryConfig = AstrometryConfig(n_det=180, n_ref=380, neighbor_k=8, eps_desc=0.015, bin_neighbor_radius=1, n_scale_bins=6, allow_reflection=True, use_mag_signature=True, use_wcs_prior=True, scale_tolerance=0.3, enforce_parity=True, stage1_n_det=50, stage1_radius_arcsec=20.0, stage2_radius_arcsec=8.0, top_k_hypotheses=60, max_quads_tested=5000, max_candidates_per_bucket=60, use_multiprobe=True, sip_degree=3, refine_clip_sigma=4.0, refine_clip_sigma_start=5.0, refine_min_match_fraction=0.5, refine_max_iter=5, refine_rematch_iters=2, refine_use_all=False, refine_n_det=None, refine_n_ref=None, refine_match_radius_arcsec=None, auto_match_resolution=True, adaptive_n_retry=2, adaptive_min_inliers=12))

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; ~2–7× more accurate than SCAMP, especially for challenging conditions, though ~4× slower.

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

Use the initial WCS to filter hypotheses by scale/parity. 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
... )