Multi-Image Camera Calibration¶

This notebook demonstrates how to calibrate a shared camera model from multiple images taken with the same optical system. By solving many images and fitting a single set of intrinsics (focal length, optical center, polynomial distortion), we can achieve much better accuracy than solving each image independently.

We use 10 TESS Full Frame Images from Camera 1, CCD 1 across different observing sectors. Since the same physical CCD is used, the camera intrinsics are shared even though the spacecraft points at different parts of the sky in each sector.

The pipeline:

1. Generate a solver database from the Hipparcos catalog
2. Load all images and extract WCS ground truth from FITS headers
3. Extract centroids from each image
4. Iteratively solve and calibrate — progressively tightening parameters across all images
5. Validate per-image results against FITS WCS
6. Visualize matched stars in each sector

Setup¶

We load images from 10 TESS sectors. TESS observes each sector for ~27 days, and the same CCD can observe different sky regions across sectors as the spacecraft repoints.

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
from IPython.display import display, HTML
import datetime
import warnings
from astropy.utils.exceptions import AstropyWarning
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio

pio.renderers.default = "notebook"

sectors = [1, 2, 3, 4, 5, 6, 13, 14, 15, 17]
fnames = [f"../../data/tess_same_ccd/sector{s:02d}_cam1_ccd1.fits" for s in sectors]

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 from IPython.display import display, HTML import datetime import warnings from astropy.utils.exceptions import AstropyWarning import plotly.graph_objects as go import plotly.express as px import plotly.io as pio pio.renderers.default = "notebook" sectors = [1, 2, 3, 4, 5, 6, 13, 14, 15, 17] fnames = [f"../../data/tess_same_ccd/sector{s:02d}_cam1_ccd1.fits" for s in sectors]

Generate the Solver Database¶

We use higher database quality parameters than the basic tutorial since multi-image calibration benefits from more verification stars and denser pattern coverage.

In [2]:

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

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

Extract WCS Ground Truth¶

We read the FITS WCS from each image header to use as ground truth for validating our calibration. We also extract the mid-exposure time and compute the Earth's barycentric velocity at that time (which can be used for stellar aberration correction, though we don't apply it here).

In [3]:

ra_dec_fits = []
frame_times = []
frame_velocities = []

for fname in fnames:
    with fits.open(fname) as hdul:
        tstart = datetime.datetime.strptime(
            hdul[0].header["DATE-OBS"], "%Y-%m-%dT%H:%M:%S.%f"
        )
        tstop = datetime.datetime.strptime(
            hdul[0].header["DATE-END"], "%Y-%m-%dT%H:%M:%S.%f"
        )
        tmid = tstart + (tstop - tstart) / 2
        frame_times.append(tmid)
        frame_velocities.append(t3rs.earth_barycentric_velocity(tmid))
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", AstropyWarning)
            wcs = WCS(hdul[1].header)
        # Center pixel in full-frame coords (accounting for 44-column overscan)
        ra_dec = wcs.pixel_to_world(1024 + 44, 1024)
        ra_dec_fits.append((ra_dec.ra.degree, ra_dec.dec.degree))

ra_dec_fits = [] frame_times = [] frame_velocities = [] for fname in fnames: with fits.open(fname) as hdul: tstart = datetime.datetime.strptime( hdul[0].header["DATE-OBS"], "%Y-%m-%dT%H:%M:%S.%f" ) tstop = datetime.datetime.strptime( hdul[0].header["DATE-END"], "%Y-%m-%dT%H:%M:%S.%f" ) tmid = tstart + (tstop - tstart) / 2 frame_times.append(tmid) frame_velocities.append(t3rs.earth_barycentric_velocity(tmid)) with warnings.catch_warnings(): warnings.simplefilter("ignore", AstropyWarning) wcs = WCS(hdul[1].header) # Center pixel in full-frame coords (accounting for 44-column overscan) ra_dec = wcs.pixel_to_world(1024 + 44, 1024) ra_dec_fits.append((ra_dec.ra.degree, ra_dec.dec.degree))

Extract Centroids¶

We extract star centroids from all 10 images using the same parameters. Each image is trimmed to the 2048×2048 science region by removing the 44-column overscan.

In [4]:

cenarr = []
for fname in fnames:
    img = fits.getdata(fname, ext=1, memmap=False).astype(np.float32)
    img = img[:2048, 44:2092]
    extraction = t3rs.extract_centroids(
        img,
        sigma_threshold=300,
        min_pixels=4,
        max_pixels=10000,
        local_bg_block_size=16,
        max_elongation=6.0,
    )
    cenarr.append(extraction.centroids)

for s, cen in zip(sectors, cenarr):
    print(f"Sector {s:2d}: {len(cen)} centroids")

cenarr = [] for fname in fnames: img = fits.getdata(fname, ext=1, memmap=False).astype(np.float32) img = img[:2048, 44:2092] extraction = t3rs.extract_centroids( img, sigma_threshold=300, min_pixels=4, max_pixels=10000, local_bg_block_size=16, max_elongation=6.0, ) cenarr.append(extraction.centroids) for s, cen in zip(sectors, cenarr): print(f"Sector {s:2d}: {len(cen)} centroids")

Sector  1: 3928 centroids
Sector  2: 3933 centroids
Sector  3: 4853 centroids
Sector  4: 5227 centroids
Sector  5: 3556 centroids
Sector  6: 2979 centroids
Sector 13: 1137 centroids
Sector 14: 965 centroids
Sector 15: 1137 centroids
Sector 17: 4068 centroids

Iterative Multi-Image Solve and Calibrate¶

This is the core of the multi-image calibration pipeline. We iterate through 4 passes, each time:

1. Solve all 10 images using the current camera model
2. Calibrate a shared camera model from all solve results simultaneously

Each pass uses a tighter match radius and higher polynomial distortion order. The shared camera model improves with each pass because:

• More images provide more star pairs to constrain the distortion fit
• Different sky regions sample different parts of the focal plane
• The alternating WCS refinement + global polynomial fit converges to a consistent solution across all images

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

In [5]:

camera_model = None

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),
]

for pass_num, (mr, ri, order, fov_err) in enumerate(pass_configs, 1):
    resarr = []
    for idx, cen in enumerate(cenarr):
        resarr.append(
            solver.solve_from_centroids(
                cen,
                fov_estimate_deg=camera_model.fov_deg if camera_model else 11.8,
                fov_max_error_deg=fov_err,
                image_shape=(2048, 2048),
                match_radius=mr,
                match_threshold=1e-5,
                refine_iterations=ri,
                camera_model=camera_model,
            )
        )

    cal = solver.calibrate_camera(
        resarr, cenarr, image_shape=(2048, 2048), order=order,
    )
    camera_model = cal.camera_model

    n_solved = sum(1 for r in resarr if r is not None)
    print(
        f"Pass {pass_num}: {n_solved}/{len(sectors)} solved, "
        f"calibration RMSE {cal.rmse_before_px:.3f} -> {cal.rmse_after_px:.3f} px, "
        f"{cal.n_inliers} inliers"
    )

camera_model = None 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), ] for pass_num, (mr, ri, order, fov_err) in enumerate(pass_configs, 1): resarr = [] for idx, cen in enumerate(cenarr): resarr.append( solver.solve_from_centroids( cen, fov_estimate_deg=camera_model.fov_deg if camera_model else 11.8, fov_max_error_deg=fov_err, image_shape=(2048, 2048), match_radius=mr, match_threshold=1e-5, refine_iterations=ri, camera_model=camera_model, ) ) cal = solver.calibrate_camera( resarr, cenarr, image_shape=(2048, 2048), order=order, ) camera_model = cal.camera_model n_solved = sum(1 for r in resarr if r is not None) print( f"Pass {pass_num}: {n_solved}/{len(sectors)} solved, " f"calibration RMSE {cal.rmse_before_px:.3f} -> {cal.rmse_after_px:.3f} px, " f"{cal.n_inliers} inliers" )

Pass 1: 10/10 solved, calibration RMSE 16.872 -> 0.522 px, 2006 inliers

Pass 2: 10/10 solved, calibration RMSE 13.910 -> 0.131 px, 2350 inliers

Pass 3: 10/10 solved, calibration RMSE 8.272 -> 0.135 px, 1996 inliers

Pass 4: 10/10 solved, calibration RMSE 9.768 -> 0.183 px, 2137 inliers

Results¶

After calibration, we compare each image's solved pointing against the FITS WCS ground truth. The "vs FITS WCS" column shows the angular separation between the solver's boresight and the WCS center pixel — this should be well under 10" for a good calibration.

In [6]:

arcsec_per_px = resarr[0].fov_deg * 3600 / 2048

rows = []
for idx, res in enumerate(resarr):
    ra_dec_model = res.pixel_to_world(0, 0)
    coord_fits = SkyCoord(
        ra_dec_fits[idx][0] * u.degree, ra_dec_fits[idx][1] * u.degree
    )
    coord_model = SkyCoord(ra_dec_model[0] * u.degree, ra_dec_model[1] * u.degree)
    separation = coord_fits.separation(coord_model).arcsecond
    rmse_px = res.rmse_arcsec / arcsec_per_px
    rows.append(
        (
            sectors[idx],
            res.num_matches,
            f"{res.rmse_arcsec:.2f}",
            f"{rmse_px:.3f}",
            f"{separation:.2f}",
            f"{res.solve_time_ms:.1f}",
        )
    )

header = '<tr><th>Sector</th><th>Matches</th><th>RMSE (")</th><th>RMSE (px)</th><th>vs FITS WCS (")</th><th>Solve (ms)</th></tr>'
body = "".join(
    f"<tr><td>{s}</td><td>{n}</td><td>{r}</td><td>{rpx}</td><td>{sep}</td><td>{t}</td></tr>"
    for s, n, r, rpx, sep, t in rows
)
display(HTML(f"<table>{header}{body}</table>"))

arcsec_per_px = resarr[0].fov_deg * 3600 / 2048 rows = [] for idx, res in enumerate(resarr): ra_dec_model = res.pixel_to_world(0, 0) coord_fits = SkyCoord( ra_dec_fits[idx][0] * u.degree, ra_dec_fits[idx][1] * u.degree ) coord_model = SkyCoord(ra_dec_model[0] * u.degree, ra_dec_model[1] * u.degree) separation = coord_fits.separation(coord_model).arcsecond rmse_px = res.rmse_arcsec / arcsec_per_px rows.append( ( sectors[idx], res.num_matches, f"{res.rmse_arcsec:.2f}", f"{rmse_px:.3f}", f"{separation:.2f}", f"{res.solve_time_ms:.1f}", ) ) header = 'SectorMatchesRMSE (")RMSE (px)vs FITS WCS (")Solve (ms)' body = "".join( f"{s}{n}{r}{rpx}{sep}{t}" for s, n, r, rpx, sep, t in rows ) display(HTML(f"{header}{body}"))

Sector	Matches	RMSE (")	RMSE (px)	vs FITS WCS (")	Solve (ms)
1	356	2.73	0.131	0.83	54.1
2	315	2.78	0.134	0.92	48.8
3	314	3.02	0.146	0.70	46.6
4	279	2.98	0.144	0.97	45.7
5	268	2.86	0.138	0.84	41.7
6	431	2.88	0.139	0.81	57.0
13	306	3.55	0.171	1.60	26.6
14	335	2.98	0.143	1.45	32.0
15	325	5.94	0.286	3.19	25.1
17	326	21.43	1.033	1.58	44.9

Visualization¶

Each panel shows one TESS sector with matched stars overlaid:

• Cyan crosses — expected star positions projected through the solved WCS
• Yellow ellipses — detected centroid covariance ellipses for matched stars

Close agreement between the crosses and ellipses confirms that the shared camera model accurately maps pixel positions to sky coordinates across all 10 sectors.

In [7]:

from plotly.subplots import make_subplots

ncols, nrows = 2, 5
ds = 4  # downsample factor
theta = np.linspace(0, 2 * np.pi, 100)

fig = make_subplots(
    rows=nrows,
    cols=ncols,
    subplot_titles=[f"Sector {s}" for s in sectors],
    vertical_spacing=0.03,
    horizontal_spacing=0.02,
)

for idx, fname in enumerate(fnames):
    row = idx // ncols + 1
    col = idx % ncols + 1

    img = fits.getdata(fname, ext=1, memmap=False)
    img = img[:2048, 44:2092].astype(np.float32)
    vmin, vmax = np.percentile(img, 1), np.percentile(img, 99.2)
    img_ds = img[::ds, ::ds]

    img_fig = px.imshow(
        img_ds,
        binary_string=True,
        color_continuous_scale="gray",
        zmin=vmin,
        zmax=vmax,
    )
    fig.add_trace(img_fig.data[0], row=row, col=col)

    res = resarr[idx]
    cen_indices = res.matched_centroids
    cat_ids = res.matched_catalog_ids

    # Project catalog stars to pixel coordinates
    star_ra = np.array([solver.get_star_by_id(int(cid)).ra_deg for cid in cat_ids])
    star_dec = np.array([solver.get_star_by_id(int(cid)).dec_deg for cid in cat_ids])
    star_xy = res.world_to_pixel(star_ra, star_dec)

    sx = (star_xy[0] + 1024) / ds
    sy = (star_xy[1] + 1024) / ds
    fig.add_trace(
        go.Scatter(
            x=sx, y=sy, mode="markers",
            marker=dict(color="cyan", symbol="cross", size=2, line=dict(width=0.25)),
            showlegend=False, hoverinfo="skip",
        ),
        row=row, col=col,
    )

    # Draw centroid covariance ellipses
    for ci in cen_indices:
        c = cenarr[idx][int(ci)]
        cx = (c.x + 1024) / ds
        cy = (c.y + 1024) / ds
        a = 6 * np.sqrt(c.cov[0, 0]) / ds
        b = 6 * np.sqrt(c.cov[1, 1]) / ds
        angle = (
            0.5 * np.arctan2(2 * c.cov[1, 0], c.cov[0, 0] - c.cov[1, 1]) * 180 / np.pi
            + 90
        )
        x = a * np.cos(theta)
        y = b * np.sin(theta)
        x_rot = x * np.cos(np.radians(angle)) - y * np.sin(np.radians(angle))
        y_rot = x * np.sin(np.radians(angle)) + y * np.cos(np.radians(angle))
        fig.add_trace(
            go.Scatter(
                x=cx + x_rot, y=cy + y_rot, mode="lines",
                line=dict(color="yellow", width=1),
                showlegend=False, hoverinfo="skip",
            ),
            row=row, col=col,
        )

fig.update_xaxes(showticklabels=False)
fig.update_yaxes(showticklabels=False)
fig.update_xaxes(range=[0, 2048 / ds], scaleanchor="y", scaleratio=1)
fig.update_yaxes(range=[2048 / ds, 0], scaleanchor="x", scaleratio=1)
fig.update_layout(
    width=700,
    height=700 * 5 / 2,
    showlegend=False,
    title_text="TESS Images with Matched Centroids",
    margin=dict(l=10, r=10, b=10),
)
fig.show()

from plotly.subplots import make_subplots ncols, nrows = 2, 5 ds = 4 # downsample factor theta = np.linspace(0, 2 * np.pi, 100) fig = make_subplots( rows=nrows, cols=ncols, subplot_titles=[f"Sector {s}" for s in sectors], vertical_spacing=0.03, horizontal_spacing=0.02, ) for idx, fname in enumerate(fnames): row = idx // ncols + 1 col = idx % ncols + 1 img = fits.getdata(fname, ext=1, memmap=False) img = img[:2048, 44:2092].astype(np.float32) vmin, vmax = np.percentile(img, 1), np.percentile(img, 99.2) img_ds = img[::ds, ::ds] img_fig = px.imshow( img_ds, binary_string=True, color_continuous_scale="gray", zmin=vmin, zmax=vmax, ) fig.add_trace(img_fig.data[0], row=row, col=col) res = resarr[idx] cen_indices = res.matched_centroids cat_ids = res.matched_catalog_ids # Project catalog stars to pixel coordinates star_ra = np.array([solver.get_star_by_id(int(cid)).ra_deg for cid in cat_ids]) star_dec = np.array([solver.get_star_by_id(int(cid)).dec_deg for cid in cat_ids]) star_xy = res.world_to_pixel(star_ra, star_dec) sx = (star_xy[0] + 1024) / ds sy = (star_xy[1] + 1024) / ds fig.add_trace( go.Scatter( x=sx, y=sy, mode="markers", marker=dict(color="cyan", symbol="cross", size=2, line=dict(width=0.25)), showlegend=False, hoverinfo="skip", ), row=row, col=col, ) # Draw centroid covariance ellipses for ci in cen_indices: c = cenarr[idx][int(ci)] cx = (c.x + 1024) / ds cy = (c.y + 1024) / ds a = 6 * np.sqrt(c.cov[0, 0]) / ds b = 6 * np.sqrt(c.cov[1, 1]) / ds angle = ( 0.5 * np.arctan2(2 * c.cov[1, 0], c.cov[0, 0] - c.cov[1, 1]) * 180 / np.pi + 90 ) x = a * np.cos(theta) y = b * np.sin(theta) x_rot = x * np.cos(np.radians(angle)) - y * np.sin(np.radians(angle)) y_rot = x * np.sin(np.radians(angle)) + y * np.cos(np.radians(angle)) fig.add_trace( go.Scatter( x=cx + x_rot, y=cy + y_rot, mode="lines", line=dict(color="yellow", width=1), showlegend=False, hoverinfo="skip", ), row=row, col=col, ) fig.update_xaxes(showticklabels=False) fig.update_yaxes(showticklabels=False) fig.update_xaxes(range=[0, 2048 / ds], scaleanchor="y", scaleratio=1) fig.update_yaxes(range=[2048 / ds, 0], scaleanchor="x", scaleratio=1) fig.update_layout( width=700, height=700 * 5 / 2, showlegend=False, title_text="TESS Images with Matched Centroids", margin=dict(l=10, r=10, b=10), ) fig.show()