{
"cells": [
{
"cell_type": "markdown",
"id": "pgv2wtx6t3",
"metadata": {},
"source": "# Basic Plate Solving\n\nThis notebook demonstrates how to use tetra3rs to solve a single star field image.\nWe'll walk through the full pipeline:\n\n1. **Generate a solver database** from the Gaia catalog\n2. **Load and prepare** a TESS Full Frame Image\n3. **Extract star centroids** from the image\n4. **Iteratively solve and calibrate** the camera model with progressively tighter parameters\n5. **Validate** the solution against the FITS WCS header\n6. **Visualize** matched stars overlaid on the image"
},
{
"cell_type": "markdown",
"id": "b1f9a001",
"metadata": {},
"source": [
"## Setup"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "73d75591",
"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",
"import plotly.graph_objects as go\n",
"import plotly.express as px\n",
"import plotly.io as pio\n",
"\n",
"pio.renderers.default = \"notebook\""
]
},
{
"cell_type": "markdown",
"id": "a2b1c3d4",
"metadata": {},
"source": "## Generate the Solver Database\n\nFirst, we build a solver database from the Gaia DR3 catalog (bundled\nautomatically via the `gaia-catalog` package — no manual download needed).\nThis precomputes geometric hash patterns for star quads, which enables fast\nplate solving. Key parameters:\n\n- `max_fov_deg=15.0` — maximum field of view the database will support\n- `pattern_max_error=0.005` — tolerance for pattern matching\n- `verification_stars_per_fov=1000` — catalog stars available per FOV for match verification\n- `epoch_proper_motion_year=2018` — propagate star positions to the TESS observation epoch"
},
{
"cell_type": "code",
"execution_count": null,
"id": "0a210956",
"metadata": {},
"outputs": [],
"source": [
"solver = t3rs.SolverDatabase.generate_from_gaia(\n",
" max_fov_deg=15.0,\n",
" pattern_max_error=0.005,\n",
" lattice_field_oversampling=100,\n",
" patterns_per_lattice_field=100,\n",
" verification_stars_per_fov=1000,\n",
" epoch_proper_motion_year=2018,\n",
")"
]
},
{
"cell_type": "markdown",
"id": "c3d4e5f6",
"metadata": {},
"source": [
"## Load and Prepare the Image\n",
"\n",
"We load a TESS Full Frame Image (Camera 1, CCD 1, Sector 1). TESS FFIs have the\n",
"science data in HDU 1 with extra columns for calibration. We trim to the\n",
"2048×2048 science region by removing the 44-column overscan on the left."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d4e5f6a7",
"metadata": {},
"outputs": [],
"source": [
"fname = \"../../data/tess_same_ccd/sector01_cam1_ccd1.fits\"\n",
"img = fits.getdata(fname, ext=1, memmap=False).astype(np.float32)\n",
"header = fits.getheader(fname, ext=1)\n",
"\n",
"img = img[:2048, 44:2092] # 2048x2048 science region\n",
"print(f\"Image shape: {img.shape}, dtype: {img.dtype}\")"
]
},
{
"cell_type": "markdown",
"id": "e5f6a7b8",
"metadata": {},
"source": [
"## Extract Star Centroids\n",
"\n",
"`extract_centroids` detects stars by finding connected pixel regions above a\n",
"local background threshold, then computes the intensity-weighted centroid\n",
"and covariance for each. Centroids are returned in **centered pixel coordinates**\n",
"(origin at image center, +X right, +Y down).\n",
"\n",
"The extraction parameters are tuned for TESS's wide-field, defocused optics:\n",
"- `sigma_threshold=300` — high threshold to reject background noise in TESS's crowded fields\n",
"- `min_pixels=4`, `max_pixels=10000` — TESS PSFs are large and defocused\n",
"- `max_elongation=6.0` — reject elongated artifacts while keeping slightly trailed stars"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f6a7b8c9",
"metadata": {},
"outputs": [],
"source": [
"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",
"centroids = extraction.centroids\n",
"print(\n",
" f\"{len(centroids)} centroids extracted (from {extraction.num_blobs_raw} raw blobs)\"\n",
")\n",
"print(\n",
" f\"Background: mean={extraction.background_mean:.1f}, sigma={extraction.background_sigma:.1f}\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "a7b8c9d0",
"metadata": {},
"source": [
"## Extract WCS Ground Truth from FITS Header\n",
"\n",
"The TESS FITS header contains a WCS solution that we'll use to validate our\n",
"plate solve. We extract the RA/Dec of the reference pixel (CRPIX) and convert\n",
"it to centered image coordinates to compare against our solver output."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b8c9d0e1",
"metadata": {},
"outputs": [],
"source": [
"import warnings\n",
"from astropy.utils.exceptions import AstropyWarning\n",
"\n",
"with warnings.catch_warnings():\n",
" warnings.simplefilter(\"ignore\", AstropyWarning)\n",
" wcs = WCS(header)\n",
"\n",
"# CRPIX in 0-indexed FITS coordinates\n",
"crpix_x = wcs.wcs.crpix[0] - 1.0\n",
"crpix_y = wcs.wcs.crpix[1] - 1.0\n",
"sky_crpix = wcs.pixel_to_world(crpix_x, crpix_y)\n",
"fits_ra = sky_crpix.ra.deg\n",
"fits_dec = sky_crpix.dec.deg\n",
"\n",
"# Convert CRPIX to centered image coordinates (accounting for the 44-column crop)\n",
"col_offset = 44\n",
"crpix_cx = (crpix_x - col_offset) - img.shape[1] / 2.0\n",
"crpix_cy = crpix_y - img.shape[0] / 2.0\n",
"\n",
"print(f\"FITS WCS reference point: RA={fits_ra:.4f}, Dec={fits_dec:.4f}\")"
]
},
{
"cell_type": "markdown",
"id": "c9d0e1f2",
"metadata": {},
"source": [
"## Iterative Solve and Calibrate\n",
"\n",
"We solve the image in multiple passes, progressively tightening the match radius\n",
"and increasing the polynomial distortion order. Each pass:\n",
"\n",
"1. **Solve** — find the camera attitude using the current camera model\n",
"2. **Calibrate** — fit a camera model (focal length, optical center, polynomial distortion)\n",
" from the matched star pairs\n",
"\n",
"The calibrated camera model from each pass feeds into the next, allowing the\n",
"solver to match more stars at tighter tolerances as the distortion model improves.\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": "d0e1f2a3",
"metadata": {},
"outputs": [],
"source": [
"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",
"camera_model = None\n",
"fov_estimate = 12.0 # approximate TESS CCD FOV in degrees\n",
"fits_coord = SkyCoord(ra=fits_ra * u.deg, dec=fits_dec * u.deg)\n",
"\n",
"for pass_num, (mr, ri, order, fov_err) in enumerate(pass_configs, 1):\n",
" result = solver.solve_from_centroids(\n",
" centroids,\n",
" fov_estimate_deg=fov_estimate,\n",
" fov_max_error_deg=fov_err,\n",
" image_shape=img.shape,\n",
" match_radius=mr,\n",
" match_threshold=1e-5,\n",
" refine_iterations=ri,\n",
" camera_model=camera_model,\n",
" )\n",
"\n",
" if result is None:\n",
" print(f\"Pass {pass_num}: FAILED to solve\")\n",
" break\n",
"\n",
" fov_estimate = result.fov_deg\n",
"\n",
" # Compare solver solution against FITS WCS at the CRPIX reference point\n",
" pred_ra, pred_dec = result.pixel_to_world(crpix_cx, crpix_cy)\n",
" pred_coord = SkyCoord(ra=pred_ra * u.deg, dec=pred_dec * u.deg)\n",
" sep = pred_coord.separation(fits_coord).to(u.arcsec)\n",
"\n",
" print(\n",
" f\"Pass {pass_num}: {result.num_matches} matches, \"\n",
" f'RMSE={result.rmse_arcsec:.2f}\", '\n",
" f\"vs FITS WCS={sep:.2f}\"\n",
" )\n",
"\n",
" cal = solver.calibrate_camera(\n",
" result,\n",
" centroids,\n",
" image_shape=img.shape,\n",
" order=order,\n",
" )\n",
" camera_model = cal.camera_model\n",
"\n",
" print(\n",
" f\" Calibration: RMSE {cal.rmse_before_px:.3f} -> {cal.rmse_after_px:.3f} px, \"\n",
" f\"{cal.n_inliers} inliers, {cal.n_outliers} outliers\"\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "e1f2a3b4",
"metadata": {},
"source": [
"## Results Summary\n",
"\n",
"After iterative calibration, the final solve result contains the camera attitude,\n",
"matched star count, and residual statistics. We compare the solved pointing\n",
"against the FITS WCS to confirm agreement."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f2a3b4c5",
"metadata": {},
"outputs": [],
"source": [
"pred_ra, pred_dec = result.pixel_to_world(crpix_cx, crpix_cy)\n",
"pred_coord = SkyCoord(ra=pred_ra * u.deg, dec=pred_dec * u.deg)\n",
"sep = pred_coord.separation(fits_coord)\n",
"\n",
"arcsec_per_px = result.fov_deg * 3600 / img.shape[1]\n",
"\n",
"print(f\"Boresight: RA={result.ra_deg:.4f}, Dec={result.dec_deg:.4f}\")\n",
"print(f\"FOV: {result.fov_deg:.4f} deg\")\n",
"print(f\"Pixel scale: {arcsec_per_px:.2f} arcsec/px\")\n",
"print(f\"Matched stars: {result.num_matches}\")\n",
"print(\n",
" f'RMSE: {result.rmse_arcsec:.2f}\" ({result.rmse_arcsec / arcsec_per_px:.3f} px)'\n",
")\n",
"print(f\"vs FITS WCS (CRPIX): {sep.to(u.arcsec):.2f}\")"
]
},
{
"cell_type": "markdown",
"id": "a3b4c5d6",
"metadata": {},
"source": [
"## Visualization\n",
"\n",
"We overlay the extracted centroids on the TESS image:\n",
"- **Yellow ellipses** — centroids matched to catalog stars\n",
"- **Red ellipses** — unmatched centroids\n",
"- **Cyan lines** — residual vectors from each matched centroid to its expected\n",
" position according to the FITS WCS (showing the agreement between solutions)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b4c5d6e7",
"metadata": {},
"outputs": [],
"source": [
"def add_ellipse(fig, cx, cy, a, b, angle_deg, **kwargs):\n",
" \"\"\"Add a rotated ellipse trace to a Plotly figure.\"\"\"\n",
" theta = np.linspace(0, 2 * np.pi, 100)\n",
" angle = np.deg2rad(angle_deg)\n",
" x = a * np.cos(theta)\n",
" y = b * np.sin(theta)\n",
" x_rot = cx + x * np.cos(angle) - y * np.sin(angle)\n",
" y_rot = cy + x * np.sin(angle) + y * np.cos(angle)\n",
" fig.add_trace(\n",
" go.Scatter(\n",
" x=x_rot,\n",
" y=y_rot,\n",
" mode=\"lines\",\n",
" showlegend=False,\n",
" hoverinfo=\"skip\",\n",
" **kwargs,\n",
" )\n",
" )\n",
"\n",
"\n",
"# Create the base image\n",
"fig = px.imshow(\n",
" img / np.max(img) * 50,\n",
" color_continuous_scale=\"gray\",\n",
" zmin=0,\n",
" zmax=1,\n",
" binary_string=True,\n",
" width=800,\n",
" height=800,\n",
" title=\"TESS Sector 1 — Matched Stars\",\n",
")\n",
"\n",
"# Draw centroid ellipses (3-sigma covariance)\n",
"matched_idx = set(int(i) for i in result.matched_centroids)\n",
"\n",
"for ix, c in enumerate(centroids):\n",
" cx = c.x + img.shape[1] / 2\n",
" cy = c.y + img.shape[0] / 2\n",
" if c.cov is not None:\n",
" a = 6 * np.sqrt(c.cov[0, 0])\n",
" b = 6 * np.sqrt(c.cov[1, 1])\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",
" else:\n",
" a, b, angle = 3, 3, 0\n",
"\n",
" is_matched = ix in matched_idx\n",
" add_ellipse(\n",
" fig,\n",
" cx,\n",
" cy,\n",
" a,\n",
" b,\n",
" angle,\n",
" line=dict(\n",
" color=\"yellow\" if is_matched else \"red\",\n",
" width=1.0 if is_matched else 0.25,\n",
" ),\n",
" )\n",
"\n",
"# Draw residual lines from matched centroids to FITS WCS expected positions\n",
"line_x, line_y = [], []\n",
"for cid, cidx in zip(result.matched_catalog_ids, result.matched_centroids):\n",
" star = solver.get_star_by_id(int(cid))\n",
" if star is None:\n",
" continue\n",
" c = centroids[int(cidx)]\n",
" act_x = c.x + img.shape[1] / 2\n",
" act_y = c.y + img.shape[0] / 2\n",
"\n",
" sky = SkyCoord(ra=star.ra_deg * u.deg, dec=star.dec_deg * u.deg)\n",
" fits_px, fits_py = wcs.world_to_pixel(sky)\n",
" exp_x = float(fits_px) - col_offset\n",
" exp_y = float(fits_py)\n",
"\n",
" line_x += [act_x, exp_x, None]\n",
" line_y += [act_y, exp_y, None]\n",
"\n",
"fig.add_trace(\n",
" go.Scatter(\n",
" x=line_x,\n",
" y=line_y,\n",
" mode=\"lines\",\n",
" line=dict(color=\"cyan\", width=1),\n",
" showlegend=False,\n",
" )\n",
")\n",
"\n",
"fig.update_xaxes(range=[0, img.shape[1]], scaleanchor=\"y\", scaleratio=1)\n",
"fig.update_yaxes(range=[img.shape[0], 0], scaleanchor=\"x\", scaleratio=1)\n",
"fig.update_coloraxes(showscale=False)\n",
"fig.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}