Basic Plate Solving¶

This notebook demonstrates how to use tetra3rs to solve a single star field image. We'll walk through the full pipeline:

1. Generate a solver database from the Hipparcos catalog
2. Load and prepare a TESS Full Frame Image
3. Extract star centroids from the image
4. Iteratively solve and calibrate the camera model with progressively tighter parameters
5. Validate the solution against the FITS WCS header
6. Visualize matched stars overlaid on the image

Setup¶

In [1]:

import tetra3rs as t3rs
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
import astropy.units as u
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio

pio.renderers.default = "notebook"

import tetra3rs as t3rs import numpy as np from astropy.io import fits from astropy.wcs import WCS from astropy.coordinates import SkyCoord import astropy.units as u import plotly.graph_objects as go import plotly.express as px import plotly.io as pio pio.renderers.default = "notebook"

Generate the Solver Database¶

First, we build a solver database from the Hipparcos-2 catalog (hip2.dat). This precomputes geometric hash patterns for star quads, which enables fast plate solving. Key parameters:

• max_fov_deg=15.0 — maximum field of view the database will support
• pattern_max_error=0.005 — tolerance for pattern matching
• verification_stars_per_fov=1000 — catalog stars available per FOV for match verification
• epoch_proper_motion_year=2018 — propagate star positions to the TESS observation epoch

In [2]:

solver = t3rs.SolverDatabase.generate_from_hipparcos(
    "../../data/hip2.dat",
    max_fov_deg=15.0,
    pattern_max_error=0.005,
    lattice_field_oversampling=100,
    patterns_per_lattice_field=100,
    verification_stars_per_fov=1000,
    epoch_proper_motion_year=2018,
)

solver = t3rs.SolverDatabase.generate_from_hipparcos( "../../data/hip2.dat", max_fov_deg=15.0, pattern_max_error=0.005, lattice_field_oversampling=100, patterns_per_lattice_field=100, verification_stars_per_fov=1000, epoch_proper_motion_year=2018, )

Load and Prepare the Image¶

We load a TESS Full Frame Image (Camera 1, CCD 1, Sector 1). TESS FFIs have the science data in HDU 1 with extra columns for calibration. We trim to the 2048×2048 science region by removing the 44-column overscan on the left.

In [3]:

fname = "../../data/tess_same_ccd/sector01_cam1_ccd1.fits"
img = fits.getdata(fname, ext=1, memmap=False).astype(np.float32)
header = fits.getheader(fname, ext=1)

img = img[:2048, 44:2092]  # 2048x2048 science region
print(f"Image shape: {img.shape}, dtype: {img.dtype}")

fname = "../../data/tess_same_ccd/sector01_cam1_ccd1.fits" img = fits.getdata(fname, ext=1, memmap=False).astype(np.float32) header = fits.getheader(fname, ext=1) img = img[:2048, 44:2092] # 2048x2048 science region print(f"Image shape: {img.shape}, dtype: {img.dtype}")

Image shape: (2048, 2048), dtype: float32

Extract Star Centroids¶

extract_centroids detects stars by finding connected pixel regions above a local background threshold, then computes the intensity-weighted centroid and covariance for each. Centroids are returned in centered pixel coordinates (origin at image center, +X right, +Y down).

The extraction parameters are tuned for TESS's wide-field, defocused optics:

• sigma_threshold=300 — high threshold to reject background noise in TESS's crowded fields
• min_pixels=4, max_pixels=10000 — TESS PSFs are large and defocused
• max_elongation=6.0 — reject elongated artifacts while keeping slightly trailed stars

In [4]:

extraction = t3rs.extract_centroids(
    img,
    sigma_threshold=300,
    min_pixels=4,
    max_pixels=10000,
    local_bg_block_size=16,
    max_elongation=6.0,
)
centroids = extraction.centroids
print(f"{len(centroids)} centroids extracted (from {extraction.num_blobs_raw} raw blobs)")
print(f"Background: mean={extraction.background_mean:.1f}, sigma={extraction.background_sigma:.1f}")

extraction = t3rs.extract_centroids( img, sigma_threshold=300, min_pixels=4, max_pixels=10000, local_bg_block_size=16, max_elongation=6.0, ) centroids = extraction.centroids print(f"{len(centroids)} centroids extracted (from {extraction.num_blobs_raw} raw blobs)") print(f"Background: mean={extraction.background_mean:.1f}, sigma={extraction.background_sigma:.1f}")

3928 centroids extracted (from 3928 raw blobs)
Background: mean=-0.4, sigma=2.0

Extract WCS Ground Truth from FITS Header¶

The TESS FITS header contains a WCS solution that we'll use to validate our plate solve. We extract the RA/Dec of the reference pixel (CRPIX) and convert it to centered image coordinates to compare against our solver output.

In [5]:

import warnings
from astropy.utils.exceptions import AstropyWarning

with warnings.catch_warnings():
    warnings.simplefilter("ignore", AstropyWarning)
    wcs = WCS(header)

# CRPIX in 0-indexed FITS coordinates
crpix_x = wcs.wcs.crpix[0] - 1.0
crpix_y = wcs.wcs.crpix[1] - 1.0
sky_crpix = wcs.pixel_to_world(crpix_x, crpix_y)
fits_ra = sky_crpix.ra.deg
fits_dec = sky_crpix.dec.deg

# Convert CRPIX to centered image coordinates (accounting for the 44-column crop)
col_offset = 44
crpix_cx = (crpix_x - col_offset) - img.shape[1] / 2.0
crpix_cy = crpix_y - img.shape[0] / 2.0

print(f"FITS WCS reference point: RA={fits_ra:.4f}, Dec={fits_dec:.4f}")

import warnings from astropy.utils.exceptions import AstropyWarning with warnings.catch_warnings(): warnings.simplefilter("ignore", AstropyWarning) wcs = WCS(header) # CRPIX in 0-indexed FITS coordinates crpix_x = wcs.wcs.crpix[0] - 1.0 crpix_y = wcs.wcs.crpix[1] - 1.0 sky_crpix = wcs.pixel_to_world(crpix_x, crpix_y) fits_ra = sky_crpix.ra.deg fits_dec = sky_crpix.dec.deg # Convert CRPIX to centered image coordinates (accounting for the 44-column crop) col_offset = 44 crpix_cx = (crpix_x - col_offset) - img.shape[1] / 2.0 crpix_cy = crpix_y - img.shape[0] / 2.0 print(f"FITS WCS reference point: RA={fits_ra:.4f}, Dec={fits_dec:.4f}")

FITS WCS reference point: RA=319.4032, Dec=-41.2818

Iterative Solve and Calibrate¶

We solve the image in multiple passes, progressively tightening the match radius and increasing the polynomial distortion order. Each pass:

1. Solve — find the camera attitude using the current camera model
2. Calibrate — fit a camera model (focal length, optical center, polynomial distortion) from the matched star pairs

The calibrated camera model from each pass feeds into the next, allowing the solver to match more stars at tighter tolerances as the distortion model improves.

Pass	Match Radius	Distortion Order
1	0.01	3
2	0.005	4
3	0.003	5
4	0.002	6

In [6]:

pass_configs = [
    # (match_radius, refine_iterations, distortion_order, fov_error_deg)
    (0.01,  10, 3, 0.5),
    (0.005, 10, 4, 0.5),
    (0.003, 10, 5, 0.5),
    (0.002, 10, 6, 0.5),
]

camera_model = None
fov_estimate = 12.0  # approximate TESS CCD FOV in degrees
fits_coord = SkyCoord(ra=fits_ra * u.deg, dec=fits_dec * u.deg)

for pass_num, (mr, ri, order, fov_err) in enumerate(pass_configs, 1):
    result = solver.solve_from_centroids(
        centroids,
        fov_estimate_deg=fov_estimate,
        fov_max_error_deg=fov_err,
        image_shape=img.shape,
        match_radius=mr,
        match_threshold=1e-5,
        refine_iterations=ri,
        camera_model=camera_model,
    )

    if result is None:
        print(f"Pass {pass_num}: FAILED to solve")
        break

    fov_estimate = result.fov_deg

    # Compare solver solution against FITS WCS at the CRPIX reference point
    pred_ra, pred_dec = result.pixel_to_world(crpix_cx, crpix_cy)
    pred_coord = SkyCoord(ra=pred_ra * u.deg, dec=pred_dec * u.deg)
    sep = pred_coord.separation(fits_coord).to(u.arcsec)

    print(
        f"Pass {pass_num}: {result.num_matches} matches, "
        f'RMSE={result.rmse_arcsec:.2f}", '
        f'vs FITS WCS={sep:.2f}'
    )

    cal = solver.calibrate_camera(
        result, centroids, image_shape=img.shape, order=order,
    )
    camera_model = cal.camera_model

    print(
        f"  Calibration: RMSE {cal.rmse_before_px:.3f} -> {cal.rmse_after_px:.3f} px, "
        f"{cal.n_inliers} inliers, {cal.n_outliers} outliers"
    )

pass_configs = [ # (match_radius, refine_iterations, distortion_order, fov_error_deg) (0.01, 10, 3, 0.5), (0.005, 10, 4, 0.5), (0.003, 10, 5, 0.5), (0.002, 10, 6, 0.5), ] camera_model = None fov_estimate = 12.0 # approximate TESS CCD FOV in degrees fits_coord = SkyCoord(ra=fits_ra * u.deg, dec=fits_dec * u.deg) for pass_num, (mr, ri, order, fov_err) in enumerate(pass_configs, 1): result = solver.solve_from_centroids( centroids, fov_estimate_deg=fov_estimate, fov_max_error_deg=fov_err, image_shape=img.shape, match_radius=mr, match_threshold=1e-5, refine_iterations=ri, camera_model=camera_model, ) if result is None: print(f"Pass {pass_num}: FAILED to solve") break fov_estimate = result.fov_deg # Compare solver solution against FITS WCS at the CRPIX reference point pred_ra, pred_dec = result.pixel_to_world(crpix_cx, crpix_cy) pred_coord = SkyCoord(ra=pred_ra * u.deg, dec=pred_dec * u.deg) sep = pred_coord.separation(fits_coord).to(u.arcsec) print( f"Pass {pass_num}: {result.num_matches} matches, " f'RMSE={result.rmse_arcsec:.2f}", ' f'vs FITS WCS={sep:.2f}' ) cal = solver.calibrate_camera( result, centroids, image_shape=img.shape, order=order, ) camera_model = cal.camera_model print( f" Calibration: RMSE {cal.rmse_before_px:.3f} -> {cal.rmse_after_px:.3f} px, " f"{cal.n_inliers} inliers, {cal.n_outliers} outliers" )

Pass 1: 307 matches, RMSE=172.62", vs FITS WCS=15.74 arcsec
  Calibration: RMSE 7.658 -> 0.470 px, 274 inliers, 33 outliers
Pass 2: 306 matches, RMSE=3.30", vs FITS WCS=3.20 arcsec
  Calibration: RMSE 12.547 -> 0.121 px, 306 inliers, 0 outliers
Pass 3: 325 matches, RMSE=2.67", vs FITS WCS=3.33 arcsec
  Calibration: RMSE 16.015 -> 0.118 px, 325 inliers, 0 outliers
Pass 4: 332 matches, RMSE=2.61", vs FITS WCS=3.36 arcsec
  Calibration: RMSE 18.209 -> 0.116 px, 332 inliers, 0 outliers

Results Summary¶

After iterative calibration, the final solve result contains the camera attitude, matched star count, and residual statistics. We compare the solved pointing against the FITS WCS to confirm agreement.

In [7]:

pred_ra, pred_dec = result.pixel_to_world(crpix_cx, crpix_cy)
pred_coord = SkyCoord(ra=pred_ra * u.deg, dec=pred_dec * u.deg)
sep = pred_coord.separation(fits_coord)

arcsec_per_px = result.fov_deg * 3600 / img.shape[1]

print(f"Boresight:          RA={result.ra_deg:.4f}, Dec={result.dec_deg:.4f}")
print(f"FOV:                {result.fov_deg:.4f} deg")
print(f"Pixel scale:        {arcsec_per_px:.2f} arcsec/px")
print(f"Matched stars:      {result.num_matches}")
print(f'RMSE:               {result.rmse_arcsec:.2f}" ({result.rmse_arcsec / arcsec_per_px:.3f} px)')
print(f'vs FITS WCS (CRPIX): {sep.to(u.arcsec):.2f}')

pred_ra, pred_dec = result.pixel_to_world(crpix_cx, crpix_cy) pred_coord = SkyCoord(ra=pred_ra * u.deg, dec=pred_dec * u.deg) sep = pred_coord.separation(fits_coord) arcsec_per_px = result.fov_deg * 3600 / img.shape[1] print(f"Boresight: RA={result.ra_deg:.4f}, Dec={result.dec_deg:.4f}") print(f"FOV: {result.fov_deg:.4f} deg") print(f"Pixel scale: {arcsec_per_px:.2f} arcsec/px") print(f"Matched stars: {result.num_matches}") print(f'RMSE: {result.rmse_arcsec:.2f}" ({result.rmse_arcsec / arcsec_per_px:.3f} px)') print(f'vs FITS WCS (CRPIX): {sep.to(u.arcsec):.2f}')

Boresight:          RA=319.5283, Dec=-41.1047
FOV:                11.8359 deg
Pixel scale:        20.81 arcsec/px
Matched stars:      332
RMSE:               2.61" (0.125 px)
vs FITS WCS (CRPIX): 3.36 arcsec

Visualization¶

We overlay the extracted centroids on the TESS image:

• Yellow ellipses — centroids matched to catalog stars
• Red ellipses — unmatched centroids
• Cyan lines — residual vectors from each matched centroid to its expected position according to the FITS WCS (showing the agreement between solutions)

In [8]:

def add_ellipse(fig, cx, cy, a, b, angle_deg, **kwargs):
    """Add a rotated ellipse trace to a Plotly figure."""
    theta = np.linspace(0, 2 * np.pi, 100)
    angle = np.deg2rad(angle_deg)
    x = a * np.cos(theta)
    y = b * np.sin(theta)
    x_rot = cx + x * np.cos(angle) - y * np.sin(angle)
    y_rot = cy + x * np.sin(angle) + y * np.cos(angle)
    fig.add_trace(go.Scatter(
        x=x_rot, y=y_rot, mode="lines",
        showlegend=False, hoverinfo="skip", **kwargs,
    ))


# Create the base image
fig = px.imshow(
    img / np.max(img) * 50,
    color_continuous_scale="gray", zmin=0, zmax=1,
    binary_string=True, width=800, height=800,
    title="TESS Sector 1 — Matched Stars",
)

# Draw centroid ellipses (3-sigma covariance)
matched_idx = set(int(i) for i in result.matched_centroids)

for ix, c in enumerate(centroids):
    cx = c.x + img.shape[1] / 2
    cy = c.y + img.shape[0] / 2
    if c.cov is not None:
        a = 6 * np.sqrt(c.cov[0, 0])
        b = 6 * np.sqrt(c.cov[1, 1])
        angle = 0.5 * np.arctan2(2 * c.cov[1, 0], c.cov[0, 0] - c.cov[1, 1]) * 180 / np.pi + 90
    else:
        a, b, angle = 3, 3, 0

    is_matched = ix in matched_idx
    add_ellipse(
        fig, cx, cy, a, b, angle,
        line=dict(
            color="yellow" if is_matched else "red",
            width=1.0 if is_matched else 0.25,
        ),
    )

# Draw residual lines from matched centroids to FITS WCS expected positions
line_x, line_y = [], []
for cid, cidx in zip(result.matched_catalog_ids, result.matched_centroids):
    star = solver.get_star_by_id(int(cid))
    if star is None:
        continue
    c = centroids[int(cidx)]
    act_x = c.x + img.shape[1] / 2
    act_y = c.y + img.shape[0] / 2

    sky = SkyCoord(ra=star.ra_deg * u.deg, dec=star.dec_deg * u.deg)
    fits_px, fits_py = wcs.world_to_pixel(sky)
    exp_x = float(fits_px) - col_offset
    exp_y = float(fits_py)

    line_x += [act_x, exp_x, None]
    line_y += [act_y, exp_y, None]

fig.add_trace(go.Scatter(
    x=line_x, y=line_y, mode="lines",
    line=dict(color="cyan", width=1), showlegend=False,
))

fig.update_xaxes(range=[0, img.shape[1]], scaleanchor="y", scaleratio=1)
fig.update_yaxes(range=[img.shape[0], 0], scaleanchor="x", scaleratio=1)
fig.update_coloraxes(showscale=False)
fig.show()

def add_ellipse(fig, cx, cy, a, b, angle_deg, **kwargs): """Add a rotated ellipse trace to a Plotly figure.""" theta = np.linspace(0, 2 * np.pi, 100) angle = np.deg2rad(angle_deg) x = a * np.cos(theta) y = b * np.sin(theta) x_rot = cx + x * np.cos(angle) - y * np.sin(angle) y_rot = cy + x * np.sin(angle) + y * np.cos(angle) fig.add_trace(go.Scatter( x=x_rot, y=y_rot, mode="lines", showlegend=False, hoverinfo="skip", **kwargs, )) # Create the base image fig = px.imshow( img / np.max(img) * 50, color_continuous_scale="gray", zmin=0, zmax=1, binary_string=True, width=800, height=800, title="TESS Sector 1 — Matched Stars", ) # Draw centroid ellipses (3-sigma covariance) matched_idx = set(int(i) for i in result.matched_centroids) for ix, c in enumerate(centroids): cx = c.x + img.shape[1] / 2 cy = c.y + img.shape[0] / 2 if c.cov is not None: a = 6 * np.sqrt(c.cov[0, 0]) b = 6 * np.sqrt(c.cov[1, 1]) angle = 0.5 * np.arctan2(2 * c.cov[1, 0], c.cov[0, 0] - c.cov[1, 1]) * 180 / np.pi + 90 else: a, b, angle = 3, 3, 0 is_matched = ix in matched_idx add_ellipse( fig, cx, cy, a, b, angle, line=dict( color="yellow" if is_matched else "red", width=1.0 if is_matched else 0.25, ), ) # Draw residual lines from matched centroids to FITS WCS expected positions line_x, line_y = [], [] for cid, cidx in zip(result.matched_catalog_ids, result.matched_centroids): star = solver.get_star_by_id(int(cid)) if star is None: continue c = centroids[int(cidx)] act_x = c.x + img.shape[1] / 2 act_y = c.y + img.shape[0] / 2 sky = SkyCoord(ra=star.ra_deg * u.deg, dec=star.dec_deg * u.deg) fits_px, fits_py = wcs.world_to_pixel(sky) exp_x = float(fits_px) - col_offset exp_y = float(fits_py) line_x += [act_x, exp_x, None] line_y += [act_y, exp_y, None] fig.add_trace(go.Scatter( x=line_x, y=line_y, mode="lines", line=dict(color="cyan", width=1), showlegend=False, )) fig.update_xaxes(range=[0, img.shape[1]], scaleanchor="y", scaleratio=1) fig.update_yaxes(range=[img.shape[0], 0], scaleanchor="x", scaleratio=1) fig.update_coloraxes(showscale=False) fig.show()