Astrometric Unwarp
Remove SIP distortion from a plate-solved image, delivering a pixel-for-pixel-correct frame with a clean linear WCS.
What this does. Takes an image with a plate solution that includes SIP (Simple Imaging Polynomial) terms and resamples it onto a distortion-free tangent-plane grid. Every output pixel is projected to the sky through the linear WCS, then back into the original frame through the full SIP WCS, and sampled there. The canvas is expanded automatically so the corners the distortion pushes outward aren't clipped. The result is a brand-new image carrying a clean linear WCS with no SIP terms — ready for mosaicking, downstream analysis, or publication.
1. When you want to unwarp
SIP terms are standard output from Astrometry.net and the SASpro Plate Solver — they encode the residual distortion that a simple TAN projection can't capture, letting the solution reach sub-pixel accuracy even for wide-field or heavily aberrated optics. SIP is a mathematical description of the distortion, not a correction of it. The pixels in your image still live where they were when the sensor read them out, so if you overlay a star catalog on the raw image, the mapping involves evaluating a polynomial at every point.
Reasons to actually bake the correction into pixels:
- Mosaicking — stitching multiple frames assumes well-behaved linear geometry. SIP terms in the input panels can create mismatches at the seams.
- Reprojection and alignment — reproject-to-common-grid routines handle linear WCS cleanly and often struggle with high-order SIP, especially when the polynomial domains disagree between frames.
- Publication figures — a clean linear WCS is more reproducible and more portable. Any FITS reader that supports basic WCS will handle the output correctly, regardless of its SIP support.
- Cosmetic cleanup — on very wide-field setups (fast astrographs, Hyperstar, camera lenses) stars near the corners are measurably distorted toward the frame edge. Unwarping actually straightens them.
- Downstream tools that assume pure linear math — photometry packages, feature detectors, machine-learning pipelines — often behave better on undistorted frames.
2. Workflow
- Plate solve first. The image must have a WCS with SIP terms. Run the SASpro Plate Solver and leave SIP enabled (the default). Without a WCS, there's nothing to invert; without SIP, there's nothing to remove.
- Open the Unwarp dialog. The active image is selected automatically. If you have multiple images open, pick the right one from the Source dropdown.
- Check the detected solution in the info readout — SIP order, pixel scale, input and output dimensions, and growth factor.
- Choose interpolation — Bicubic (default, best quality), Bilinear (faster), or Nearest (for mask-like data where interpolation would be wrong).
- Decide on canvas expansion — leave on by default. Turn off only if you absolutely need the output to be the same size as the input and accept that distortion-corrected corners may fall outside the frame.
- Pick fill behavior — black (default) or NaN (preserves "no data" information through downstream processing).
- Deliver. Apply (in place) to overwrite the current doc with undo, New View for a side-by-side comparison, or Save FITS to write straight to disk.
3. The info readout
Before you commit, the dialog shows what the unwarp will do:
SIP order N · s″/px — N is the polynomial order (2 for small aberrations, 3-4 for strong wide-field distortion). s is the mean pixel scale from the linear part of the WCS.
input W × H — your source image dimensions.
output W' × H' (canvas grows ×G) — the output size after unwarping. G is the area growth factor. If expansion is off, this stays the same as input and the report says "canvas unchanged".
⚠ large expansion — if G exceeds 2.5, the dialog warns you. In practice this means the SIP polynomial is predicting very large edge distortion. It can be legitimate for extreme optics but it's also a classic sign of a bad plate solution — a polynomial that fit noise and extrapolated wildly beyond the image bounds. Worth a sanity check before you commit.
After the unwarp completes, the status line reports the measured correction:
Peak edge-to-edge correction — the maximum pixel shift observed across the border + interior sample grid, with the arcsec equivalent shown when pixel scale is known.
Median — the typical correction. For well-behaved optics this is often small (sub-pixel to a few pixels) while the peak is at the frame corners.
Final line — a reminder that re-solving the output should produce much smaller residuals (ideally a linear-only solution with no SIP needed).
4. Delivery modes
| Mode | What happens | Best for |
|---|
| Apply (in place) | The source image is replaced with the unwarped version. The source document's WCS and header are updated. Fully undoable. | Normal workflow — commit the correction and keep processing. |
| New View | Opens the unwarped image as a new document, leaving the source untouched. | Comparing before/after, keeping the SIP-tagged version as a reference, or branching your workflow. |
| Save FITS… | Writes straight to disk as a FITS file with the rewritten linear WCS. The current doc is not modified. | Exporting to another program (DS9, Aladin, AstroImageJ) or archiving a corrected frame separately. |
5. The preset drag grip
Bottom-left of the dialog is the Star-Stretch-style preset drag handle. The current settings (interpolation, expand toggle, fill behavior) form a canonical three-field preset:
{ "expand": True|False, "order": 0|1|3, "fill_nan": True|False }
Drag this handle onto:
- the canvas — creates a shortcut icon that remembers these exact settings. Double-click it later to reopen the dialog with the same values pre-filled.
- a plate-solved image directly — runs the unwarp headlessly with these settings, no dialog, result delivered in place with undo.
The headless path respects the same guards as the dialog: if the target has no WCS or no SIP terms, the drop is a logged no-op rather than an error — nothing is modified.
6. Workflow recipes
Standard wide-field pipeline. Plate solve with SIP (default order 3). Unwarp with Bicubic + Expand + black fill. Apply in place. Continue to stretching / stacking / further processing on the clean linear frame. One drag of the preset grip onto the canvas lets you repeat the exact same unwarp settings on your next frame with a double-click.
Mosaic preparation. For each panel: plate solve → unwarp → save FITS (or push through your mosaic tool directly). All panels now share a common distortion-free coordinate system, so your mosaic blender only has to handle the panel-to-panel geometry. Use NaN fill if your mosaic tool handles no-data regions.
Sanity check a suspect plate solution. If the info readout shows "⚠ large expansion" (growth factor > 2.5), don't just click Apply. Re-run the Plate Solver with a lower SIP order (2 or even TAN-only), or with a different star catalog, or on a cropped region, and compare. A good solution typically produces growth under 1.15 for astrograph-class optics, under 1.5 for fast widefield, and only exceeds 2 on genuinely extreme systems.
Comparing before and after. Use New View instead of Apply. The source document keeps its SIP terms, the new view carries the clean linear WCS. Overlay a star catalog on each and visually compare how tightly the catalog positions sit on top of the stars in the corners of the frame. Catalog overlays using the unwarped image should sit nearly dead-center on every star across the entire field.
Preserve mask data (nearest-neighbor). If you're unwarping a label mask, segmentation output, or any image whose pixel values are categorical rather than continuous (e.g. 0/1/2 representing different star types or exclusion zones), set interpolation to Nearest. Bilinear or Bicubic would blend category values and produce meaningless intermediate numbers.
7. Interpolation choice
| Method | Order | When to use |
|---|
| Bicubic (best) | 3 | Default. Produces the sharpest, cleanest result for normal astronomical imagery. Preserves stellar PSFs without ringing on typical data. |
| Bilinear | 1 | Faster, slightly softer. Good for very large images or quick previews. Still vastly better than Nearest for continuous data. |
| Nearest | 0 | Only for categorical / label data (masks, segmentation). Never produces intermediate values, so categorical meaning is preserved. |
8. Technical details
The SIP polynomial, the two-step projection pipeline, the canvas-expansion math, the band-parallel resampler, and the header-rewrite rules.
8.1 SIP polynomial — the distortion model
SIP (Shupe et al. 2005) extends a standard FITS WCS with two polynomial corrections applied in pixel space before the linear transformation runs. Given a reference pixel (CRPIX1, CRPIX2), define pixel offsets:
u = x - \mathrm{CRPIX1}, \quad v = y - \mathrm{CRPIX2}
SIP introduces two polynomials of order n (the SIP order shown in the info readout):
f(u, v) = \sum_{p+q \leq n} A_{p,q} \cdot u^p v^q
g(u, v) = \sum_{p+q \leq n} B_{p,q} \cdot u^p v^q
The corrected (distortion-removed-in-pixel-space) position is:
x_{\mathrm{corr}} = u + f(u, v) + \mathrm{CRPIX1}
y_{\mathrm{corr}} = v + g(u, v) + \mathrm{CRPIX2}
The A and B coefficients are stored in the FITS header as A_i_j and B_i_j keywords, with A_ORDER and B_ORDER giving the maximum order. For the inverse mapping (from the corrected coordinates back to raw pixel coordinates) there are optional AP_i_j and BP_i_j polynomials, which Astrometry.net typically fits alongside the forward ones. SIP order is reported as the max of A and B order:
n = \max(A_{\mathrm{order}}, B_{\mathrm{order}})
8.2 The two-step projection pipeline
The unwarp works by, for each pixel of the output image, finding where its sky position lives in the source image. Two WCS objects are involved:
- Full SIP WCS — the solution as plate-solved, with polynomial distortion terms.
- Linear WCS — the same solution with
sip = None and -SIP stripped from the CTYPEs. This is a pure tangent-plane projection with CRPIX, CRVAL, and the CD matrix (or PC+CDELT) only.
For every output pixel (x_out, y_out):
1. Project through the linear WCS forward to get sky coordinates:
(x_{\mathrm{out}}, y_{\mathrm{out}}) \to (\alpha, \delta) \; \text{via linear WCS}
2. Project back through the full SIP WCS inverse to get source pixel coordinates:
(\alpha, \delta) \to (x_{\mathrm{src}}, y_{\mathrm{src}}) \; \text{via full SIP WCS}
3. Sample the source image at (x_src, y_src) using the chosen interpolation order:
I_{\mathrm{out}}(x_o, y_o) = I_{\mathrm{src}}(x_s, y_s)
The result is an image where sampling at (x_out, y_out) and projecting through the output's linear WCS gives the same sky position as sampling at (x_src, y_src) and projecting through the full SIP WCS. Both WCSes agree about the sky — just with and without distortion in pixel space.
8.3 Canvas expansion
Removing distortion generally changes the shape of the image's footprint on the sky. The corners of a strongly distorted frame get pushed outward when the distortion is removed, so if the output canvas is kept the same size as the input, the corrected corners fall off the frame.
To find the correct canvas size, the tool samples points densely around the input image border (96 samples per side) plus a coarse 24×24 interior grid, projects each through SIP forward and linear inverse, and takes the bounding box of all the resulting positions. The output canvas is sized to fit this bounding box, and CRPIX is adjusted so the new canvas is correctly aligned to the sky:
\mathrm{CRPIX}_{\mathrm{new}} = \mathrm{CRPIX}_{\mathrm{old}} - (x_0, y_0)
where (x_0, y_0) is the lower-left offset of the bounding box relative to the original frame. The growth factor reported in the info readout is:
\mathrm{growth} = \frac{W_{\mathrm{out}} \cdot H_{\mathrm{out}}}{W_{\mathrm{in}} \cdot H_{\mathrm{in}}}
If expansion is disabled, the canvas keeps the original dimensions and (x_0, y_0) = 0, which may clip the corrected corners.
8.4 Shift measurement
After computing the output WCS, the tool measures how far each sample point moved from its source position to its new linear-WCS position. For each sample pixel (p_x, p_y) in the source:
- Project to true sky via the full SIP WCS: (α, δ) = WCS_full(p_x, p_y).
- Project back into the output's linear WCS: (q_x, q_y) = LinearWCS-1(α, δ), accounting for the origin shift.
- Compute the shift magnitude:
\Delta = \sqrt{(x_{\mathrm{linear}} - x_{\mathrm{source}})^2 + (y_{\mathrm{linear}} - y_{\mathrm{source}})^2}
Then take the median across all samples (typical correction) and the maximum (peak correction, usually at the corners). The arcsec equivalent uses the plate scale:
\Delta_{\mathrm{arcsec}} = \Delta_{\mathrm{px}} \cdot s_{\mathrm{plate}}
where splate is the mean pixel scale derived from the determinant of the CD matrix:
s_{\mathrm{plate}} = \sqrt{|\det(\mathrm{CD})|} \cdot 3600 \; \mathrm{arcsec/px}
8.5 Band-parallel resampling
The output image is divided into horizontal bands (16-256 rows depending on frame height) and the resample is dispatched across a thread pool with worker count equal to the CPU core count. Each band builds a (y, x) meshgrid for its rows, projects the full grid through linear-forward + SIP-inverse to get source coordinates, and calls scipy.ndimage.map_coordinates once per output channel with the chosen interpolation order.
WCS objects are not thread-safe. The resampler keeps a thread-local pair of deepcopied WCS objects (one linear-out, one full-SIP-in) per worker thread, created lazily on first use. Sharing a single WCS instance across threads caused silent corruption in earlier builds — the astropy WCS internals cache intermediate state during each projection call, and concurrent calls stomp on each other's cache without raising an error. The thread-local pattern avoids this at the cost of a modest memory overhead (two small Python objects per worker).
Progress is reported by rows completed rather than bands dispatched, so the progress bar is monotonic even when bands finish out of order.
8.6 Header rewrite
The output's FITS header is built by starting from the source header and stripping two categories of keyword:
- SIP keys: anything matching
A_i_j, B_i_j, AP_i_j, BP_i_j, A_ORDER, B_ORDER, AP_ORDER, BP_ORDER, or *_DMAX.
- WCS keys:
WCSAXES, CTYPE*, CUNIT*, CRPIX*, CRVAL*, CDELT*, CROTA*, CD*_*, PC*_*, PV*_*, LONPOLE, LATPOLE.
Then the new linear WCS is written in via w_out.to_header(relax=True), which produces the full standard WCS keyword set. NAXIS1 and NAXIS2 are updated to the expanded canvas dimensions. A sentinel SASUNWRP = True is added so downstream tools can detect that this image has already been unwarped.
All other header keywords — instrument name, exposure, filter, object, date, observation metadata — are passed through unchanged. The intent is that the output is a drop-in replacement for the input in any pipeline that cares about these, with only the geometry corrected.
8.7 Verification
After an unwarp, re-running the Plate Solver on the output should produce a solution with:
- Vanishingly small SIP coefficients (or SIP disabled altogether, with a pure TAN solution converging on the same CD matrix).
- Residuals (the catalog-to-image match RMS) at least as good as the original solution.
- A CD matrix identical to the output's linear WCS CD matrix.
This round-trip test is the ultimate check that the unwarp preserved the astrometry correctly:
\mathrm{residual} = \|\mathrm{WCS}^{-1}(\mathrm{WCS}(x, y)) - (x, y)\|
should be at the level of plate-solve noise (fractions of a pixel), not systematic.
Credits. SASpro Astrometric Unwarp — the two-step linear/SIP projection pipeline, the border-plus-grid canvas-sizing algorithm, the shift-measurement reporting, the thread-local-WCS band-parallel resampler, the three-mode delivery (Apply in place, New View, Save FITS), and the preset-drag-grip headless-apply integration — is original work by Franklin Marek for SetiAstro Suite Pro. Standard components credited to the community: SIP polynomial convention (Shupe, Moshir, Li, Makovoz, Narron, Hook 2005, "The SIP Convention for Representing Distortion in FITS Image Headers"); astropy for the WCS implementation (wcs.WCS with relax=True for SIP support, all_pix2world / wcs_world2pix / all_world2pix for the forward and inverse projections, wcs.utils.proj_plane_pixel_scales for pixel scale extraction); scipy for the band resampling (ndimage.map_coordinates with configurable spline order); the FITS WCS standard (Greisen & Calabretta 2002, Calabretta & Greisen 2002); and Astrometry.net (Lang et al. 2010) as the canonical producer of the SIP solutions this tool consumes.