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 Gaia 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_gaia(
    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,
    multiscale_step=1.5,
    catalog_nside=8,
)

solver = t3rs.SolverDatabase.generate_from_gaia( 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, multiscale_step=1.5, catalog_nside=8, )

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)
        # Geometric center of the trimmed 2048x2048 science region, i.e.
        # (W-1)/2 = 1023.5 — the solver's centered-coordinate origin —
        # in full-frame coords (accounting for the 44-column overscan)
        ra_dec = wcs.pixel_to_world(1023.5 + 44, 1023.5)
        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) # Geometric center of the trimmed 2048x2048 science region, i.e. # (W-1)/2 = 1023.5 — the solver's centered-coordinate origin — # in full-frame coords (accounting for the 44-column overscan) ra_dec = wcs.pixel_to_world(1023.5 + 44, 1023.5) 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=180,
        min_pixels=3,
        max_pixels=10000,
        local_bg_block_size=16,
        max_elongation=6.0,
        matched_filter_sigma=None,
    )
    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=180, min_pixels=3, max_pixels=10000, local_bg_block_size=16, max_elongation=6.0, matched_filter_sigma=None, ) cenarr.append(extraction.centroids) for s, cen in zip(sectors, cenarr): print(f"Sector {s:2d}: {len(cen)} centroids")

Sector  1: 4057 centroids
Sector  2: 4243 centroids
Sector  3: 5163 centroids
Sector  4: 5470 centroids
Sector  5: 3905 centroids
Sector  6: 3310 centroids
Sector 13: 1286 centroids
Sector 14: 1145 centroids
Sector 15: 1275 centroids
Sector 17: 4221 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

In [5]:

camera_model = None

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

for pass_num, (mr, 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,
                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)
    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, distortion_order, fov_error_deg) (0.005, 3, 0.5), (0.005, 4, 0.5), (0.003, 5, 0.5), (0.002, 6, 0.5), ] for pass_num, (mr, 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, 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) 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: 9/10 solved, calibration RMSE 16.345 -> 0.508 px, 4195 inliers

Pass 2: 10/10 solved, calibration RMSE 13.489 -> 0.089 px, 4205 inliers

Pass 3: 10/10 solved, calibration RMSE 15.252 -> 0.079 px, 4237 inliers

Pass 4: 10/10 solved, calibration RMSE 9.604 -> 0.078 px, 3772 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 = '<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>"))

Sector	Matches	RMSE (")	RMSE (px)	vs FITS WCS (")	Solve (ms)
1	916	1.03	0.050	0.15	69.1
2	632	1.07	0.052	0.51	49.7
3	542	1.09	0.053	0.31	45.7
4	618	1.42	0.069	0.23	51.2
5	866	1.08	0.052	0.36	64.0
6	1739	1.79	0.087	0.32	114.2
13	1228	4.75	0.230	0.52	70.4
14	1089	2.67	0.129	0.38	82.9
15	1209	2.68	0.130	0.52	73.6
17	1305	1.62	0.079	0.34	89.1

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
c0 = 1023.5  # centered-origin offset: geometric center (W-1)/2 of 2048
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] + c0) / ds
    sy = (star_xy[1] + c0) / 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 + c0) / ds
        cy = (c.y + c0) / 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 c0 = 1023.5 # centered-origin offset: geometric center (W-1)/2 of 2048 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] + c0) / ds sy = (star_xy[1] + c0) / 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 + c0) / ds cy = (c.y + c0) / 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()