{
"cells": [
{
"cell_type": "markdown",
"id": "0szoz12vx6mg",
"metadata": {},
"source": [
"# Multi-Image Camera Calibration\n",
"\n",
"This notebook demonstrates how to calibrate a shared camera model from multiple\n",
"images taken with the same optical system. By solving many images and fitting a\n",
"single set of intrinsics (focal length, optical center, polynomial distortion),\n",
"we can achieve much better accuracy than solving each image independently.\n",
"\n",
"We use 10 TESS Full Frame Images from Camera 1, CCD 1 across different observing\n",
"sectors. Since the same physical CCD is used, the camera intrinsics are shared\n",
"even though the spacecraft points at different parts of the sky in each sector.\n",
"\n",
"The pipeline:\n",
"\n",
"1. **Generate a solver database** from the Hipparcos catalog\n",
"2. **Load all images** and extract WCS ground truth from FITS headers\n",
"3. **Extract centroids** from each image\n",
"4. **Iteratively solve and calibrate** — progressively tightening parameters across all images\n",
"5. **Validate** per-image results against FITS WCS\n",
"6. **Visualize** matched stars in each sector"
]
},
{
"cell_type": "markdown",
"id": "a1b2c3d4",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"We load images from 10 TESS sectors. TESS observes each sector for ~27 days,\n",
"and the same CCD can observe different sky regions across sectors as the\n",
"spacecraft repoints."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "021a1b31",
"metadata": {},
"outputs": [],
"source": [
"import tetra3rs as t3rs\n",
"import numpy as np\n",
"from astropy.io import fits\n",
"from astropy.wcs import WCS\n",
"from astropy.coordinates import SkyCoord\n",
"import astropy.units as u\n",
"from IPython.display import display, HTML\n",
"import datetime\n",
"import warnings\n",
"from astropy.utils.exceptions import AstropyWarning\n",
"import plotly.graph_objects as go\n",
"import plotly.express as px\n",
"import plotly.io as pio\n",
"\n",
"pio.renderers.default = \"notebook\"\n",
"\n",
"sectors = [1, 2, 3, 4, 5, 6, 13, 14, 15, 17]\n",
"fnames = [f\"../../data/tess_same_ccd/sector{s:02d}_cam1_ccd1.fits\" for s in sectors]"
]
},
{
"cell_type": "markdown",
"id": "b2c3d4e5",
"metadata": {},
"source": [
"## Generate the Solver Database\n",
"\n",
"We use higher database quality parameters than the basic tutorial since\n",
"multi-image calibration benefits from more verification stars and denser\n",
"pattern coverage."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ebcc9f4f",
"metadata": {},
"outputs": [],
"source": "solver = t3rs.SolverDatabase.generate_from_hipparcos(\n max_fov_deg=14.0,\n pattern_max_error=0.005,\n lattice_field_oversampling=100,\n patterns_per_lattice_field=500,\n verification_stars_per_fov=3000,\n epoch_proper_motion_year=2018,\n)"
},
{
"cell_type": "markdown",
"id": "c3d4e5f6",
"metadata": {},
"source": [
"## Extract WCS Ground Truth\n",
"\n",
"We read the FITS WCS from each image header to use as ground truth for\n",
"validating our calibration. We also extract the mid-exposure time and\n",
"compute the Earth's barycentric velocity at that time (which can be used\n",
"for stellar aberration correction, though we don't apply it here)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b0c557d5",
"metadata": {},
"outputs": [],
"source": [
"ra_dec_fits = []\n",
"frame_times = []\n",
"frame_velocities = []\n",
"\n",
"for fname in fnames:\n",
" with fits.open(fname) as hdul:\n",
" tstart = datetime.datetime.strptime(\n",
" hdul[0].header[\"DATE-OBS\"], \"%Y-%m-%dT%H:%M:%S.%f\"\n",
" )\n",
" tstop = datetime.datetime.strptime(\n",
" hdul[0].header[\"DATE-END\"], \"%Y-%m-%dT%H:%M:%S.%f\"\n",
" )\n",
" tmid = tstart + (tstop - tstart) / 2\n",
" frame_times.append(tmid)\n",
" frame_velocities.append(t3rs.earth_barycentric_velocity(tmid))\n",
" with warnings.catch_warnings():\n",
" warnings.simplefilter(\"ignore\", AstropyWarning)\n",
" wcs = WCS(hdul[1].header)\n",
" # Center pixel in full-frame coords (accounting for 44-column overscan)\n",
" ra_dec = wcs.pixel_to_world(1024 + 44, 1024)\n",
" ra_dec_fits.append((ra_dec.ra.degree, ra_dec.dec.degree))"
]
},
{
"cell_type": "markdown",
"id": "d4e5f6a7",
"metadata": {},
"source": [
"## Extract Centroids\n",
"\n",
"We extract star centroids from all 10 images using the same parameters.\n",
"Each image is trimmed to the 2048×2048 science region by removing the\n",
"44-column overscan."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6fffd2a6",
"metadata": {},
"outputs": [],
"source": [
"cenarr = []\n",
"for fname in fnames:\n",
" img = fits.getdata(fname, ext=1, memmap=False).astype(np.float32)\n",
" img = img[:2048, 44:2092]\n",
" extraction = t3rs.extract_centroids(\n",
" img,\n",
" sigma_threshold=300,\n",
" min_pixels=4,\n",
" max_pixels=10000,\n",
" local_bg_block_size=16,\n",
" max_elongation=6.0,\n",
" )\n",
" cenarr.append(extraction.centroids)\n",
"\n",
"for s, cen in zip(sectors, cenarr):\n",
" print(f\"Sector {s:2d}: {len(cen)} centroids\")"
]
},
{
"cell_type": "markdown",
"id": "e5f6a7b8",
"metadata": {},
"source": [
"## Iterative Multi-Image Solve and Calibrate\n",
"\n",
"This is the core of the multi-image calibration pipeline. We iterate through\n",
"4 passes, each time:\n",
"\n",
"1. **Solve** all 10 images using the current camera model\n",
"2. **Calibrate** a shared camera model from all solve results simultaneously\n",
"\n",
"Each pass uses a tighter match radius and higher polynomial distortion order.\n",
"The shared camera model improves with each pass because:\n",
"\n",
"- More images provide more star pairs to constrain the distortion fit\n",
"- Different sky regions sample different parts of the focal plane\n",
"- The alternating WCS refinement + global polynomial fit converges to\n",
" a consistent solution across all images\n",
"\n",
"| Pass | Match Radius | Distortion Order |\n",
"|------|-------------|------------------|\n",
"| 1 | 0.01 | 3 |\n",
"| 2 | 0.005 | 4 |\n",
"| 3 | 0.003 | 5 |\n",
"| 4 | 0.002 | 6 |"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "43bcb824",
"metadata": {},
"outputs": [],
"source": [
"camera_model = None\n",
"\n",
"pass_configs = [\n",
" # (match_radius, refine_iterations, distortion_order, fov_error_deg)\n",
" (0.01, 10, 3, 0.5),\n",
" (0.005, 10, 4, 0.5),\n",
" (0.003, 10, 5, 0.5),\n",
" (0.002, 10, 6, 0.5),\n",
"]\n",
"\n",
"for pass_num, (mr, ri, order, fov_err) in enumerate(pass_configs, 1):\n",
" resarr = []\n",
" for idx, cen in enumerate(cenarr):\n",
" resarr.append(\n",
" solver.solve_from_centroids(\n",
" cen,\n",
" fov_estimate_deg=camera_model.fov_deg if camera_model else 11.8,\n",
" fov_max_error_deg=fov_err,\n",
" image_shape=(2048, 2048),\n",
" match_radius=mr,\n",
" match_threshold=1e-5,\n",
" refine_iterations=ri,\n",
" camera_model=camera_model,\n",
" )\n",
" )\n",
"\n",
" cal = solver.calibrate_camera(\n",
" resarr, cenarr, image_shape=(2048, 2048), order=order,\n",
" )\n",
" camera_model = cal.camera_model\n",
"\n",
" n_solved = sum(1 for r in resarr if r is not None)\n",
" print(\n",
" f\"Pass {pass_num}: {n_solved}/{len(sectors)} solved, \"\n",
" f\"calibration RMSE {cal.rmse_before_px:.3f} -> {cal.rmse_after_px:.3f} px, \"\n",
" f\"{cal.n_inliers} inliers\"\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "f6a7b8c9",
"metadata": {},
"source": [
"## Results\n",
"\n",
"After calibration, we compare each image's solved pointing against the\n",
"FITS WCS ground truth. The \"vs FITS WCS\" column shows the angular separation\n",
"between the solver's boresight and the WCS center pixel — this should be\n",
"well under 10\" for a good calibration."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a7b8c9d0",
"metadata": {},
"outputs": [],
"source": [
"arcsec_per_px = resarr[0].fov_deg * 3600 / 2048\n",
"\n",
"rows = []\n",
"for idx, res in enumerate(resarr):\n",
" ra_dec_model = res.pixel_to_world(0, 0)\n",
" coord_fits = SkyCoord(\n",
" ra_dec_fits[idx][0] * u.degree, ra_dec_fits[idx][1] * u.degree\n",
" )\n",
" coord_model = SkyCoord(ra_dec_model[0] * u.degree, ra_dec_model[1] * u.degree)\n",
" separation = coord_fits.separation(coord_model).arcsecond\n",
" rmse_px = res.rmse_arcsec / arcsec_per_px\n",
" rows.append(\n",
" (\n",
" sectors[idx],\n",
" res.num_matches,\n",
" f\"{res.rmse_arcsec:.2f}\",\n",
" f\"{rmse_px:.3f}\",\n",
" f\"{separation:.2f}\",\n",
" f\"{res.solve_time_ms:.1f}\",\n",
" )\n",
" )\n",
"\n",
"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>'\n",
"body = \"\".join(\n",
" f\"<tr><td>{s}</td><td>{n}</td><td>{r}</td><td>{rpx}</td><td>{sep}</td><td>{t}</td></tr>\"\n",
" for s, n, r, rpx, sep, t in rows\n",
")\n",
"display(HTML(f\"<table>{header}{body}</table>\"))"
]
},
{
"cell_type": "markdown",
"id": "b8c9d0e1",
"metadata": {},
"source": [
"## Visualization\n",
"\n",
"Each panel shows one TESS sector with matched stars overlaid:\n",
"- **Cyan crosses** — expected star positions projected through the solved WCS\n",
"- **Yellow ellipses** — detected centroid covariance ellipses for matched stars\n",
"\n",
"Close agreement between the crosses and ellipses confirms that the shared\n",
"camera model accurately maps pixel positions to sky coordinates across\n",
"all 10 sectors."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "35d1e520",
"metadata": {},
"outputs": [],
"source": [
"from plotly.subplots import make_subplots\n",
"\n",
"ncols, nrows = 2, 5\n",
"ds = 4 # downsample factor\n",
"theta = np.linspace(0, 2 * np.pi, 100)\n",
"\n",
"fig = make_subplots(\n",
" rows=nrows,\n",
" cols=ncols,\n",
" subplot_titles=[f\"Sector {s}\" for s in sectors],\n",
" vertical_spacing=0.03,\n",
" horizontal_spacing=0.02,\n",
")\n",
"\n",
"for idx, fname in enumerate(fnames):\n",
" row = idx // ncols + 1\n",
" col = idx % ncols + 1\n",
"\n",
" img = fits.getdata(fname, ext=1, memmap=False)\n",
" img = img[:2048, 44:2092].astype(np.float32)\n",
" vmin, vmax = np.percentile(img, 1), np.percentile(img, 99.2)\n",
" img_ds = img[::ds, ::ds]\n",
"\n",
" img_fig = px.imshow(\n",
" img_ds,\n",
" binary_string=True,\n",
" color_continuous_scale=\"gray\",\n",
" zmin=vmin,\n",
" zmax=vmax,\n",
" )\n",
" fig.add_trace(img_fig.data[0], row=row, col=col)\n",
"\n",
" res = resarr[idx]\n",
" cen_indices = res.matched_centroids\n",
" cat_ids = res.matched_catalog_ids\n",
"\n",
" # Project catalog stars to pixel coordinates\n",
" star_ra = np.array([solver.get_star_by_id(int(cid)).ra_deg for cid in cat_ids])\n",
" star_dec = np.array([solver.get_star_by_id(int(cid)).dec_deg for cid in cat_ids])\n",
" star_xy = res.world_to_pixel(star_ra, star_dec)\n",
"\n",
" sx = (star_xy[0] + 1024) / ds\n",
" sy = (star_xy[1] + 1024) / ds\n",
" fig.add_trace(\n",
" go.Scatter(\n",
" x=sx, y=sy, mode=\"markers\",\n",
" marker=dict(color=\"cyan\", symbol=\"cross\", size=2, line=dict(width=0.25)),\n",
" showlegend=False, hoverinfo=\"skip\",\n",
" ),\n",
" row=row, col=col,\n",
" )\n",
"\n",
" # Draw centroid covariance ellipses\n",
" for ci in cen_indices:\n",
" c = cenarr[idx][int(ci)]\n",
" cx = (c.x + 1024) / ds\n",
" cy = (c.y + 1024) / ds\n",
" a = 6 * np.sqrt(c.cov[0, 0]) / ds\n",
" b = 6 * np.sqrt(c.cov[1, 1]) / ds\n",
" angle = (\n",
" 0.5 * np.arctan2(2 * c.cov[1, 0], c.cov[0, 0] - c.cov[1, 1]) * 180 / np.pi\n",
" + 90\n",
" )\n",
" x = a * np.cos(theta)\n",
" y = b * np.sin(theta)\n",
" x_rot = x * np.cos(np.radians(angle)) - y * np.sin(np.radians(angle))\n",
" y_rot = x * np.sin(np.radians(angle)) + y * np.cos(np.radians(angle))\n",
" fig.add_trace(\n",
" go.Scatter(\n",
" x=cx + x_rot, y=cy + y_rot, mode=\"lines\",\n",
" line=dict(color=\"yellow\", width=1),\n",
" showlegend=False, hoverinfo=\"skip\",\n",
" ),\n",
" row=row, col=col,\n",
" )\n",
"\n",
"fig.update_xaxes(showticklabels=False)\n",
"fig.update_yaxes(showticklabels=False)\n",
"fig.update_xaxes(range=[0, 2048 / ds], scaleanchor=\"y\", scaleratio=1)\n",
"fig.update_yaxes(range=[2048 / ds, 0], scaleanchor=\"x\", scaleratio=1)\n",
"fig.update_layout(\n",
" width=700,\n",
" height=700 * 5 / 2,\n",
" showlegend=False,\n",
" title_text=\"TESS Images with Matched Centroids\",\n",
" margin=dict(l=10, r=10, b=10),\n",
")\n",
"fig.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.12.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}