# Reverse polynomial model for reverse astrometry

I have a polynomial model that can transform pixel X,Y coordinates to standard $$\xi, \eta$$. This model looks like this:

$$\xi/\eta = a_1 + a_2x + a_3y + a_4x^2 + a_5xy + a_6y^2 + a_7x^3 + a_8x^2y+a_9xy^2 + a_{10}y^3 + a_{11}x^4 + a_{12}x^3y + a_{13}x^2*y^2 + a_{14}xy^3 + a_{15}y^4$$

where I know every $$a_i$$ cooefficent for both $$\xi$$ and $$\eta$$. I also know equatorial coordinates of center pixel so I can use those formulas to get equatorial coordinates of given X, Y pixel:

$$\tan\Delta\alpha = \frac{\xi}{cos\delta_0 - \eta\sin\delta_0} \\ \sin\delta = \frac{\sin\delta_0 + \eta\cos\delta_0}{\sqrt{1 + \xi^2 + \eta^2}}$$

So right now I can get equatorial coordinates of given pixel. My question is how to reverse this polynomial model to get pixel X, Y coordinates from standard coordinates?

EDIT: To answer concerns about relevancy for astronomy and inversing this using CASs: this equation is somewhat standard equation for plate solving, especially for wide field of views where you can tweak $$a_i$$ to reduce distortion. I don't want to use straight mathematical approach for inversing this function, because the results are quite complicated and I feel like, with how commonly SIP and WCS are used, there should be an easier way of calculating pixel X,Y coordinates from equatorial (or standard).

EDIT2: Sorry for misunderstanding. I don't have corresponding image or valid FITS header. Only thing I have got is SIP (the $$a_i$$ coefficetns and $$\alpha_0$$ and $$\delta_0$$) model for converting from pixel to standard.

• @planetmaker any question whose answer directly helps an astrometer perform astrometric calculations is 100% relevant to astronomy and is 100% on topic. That's prima facie and tautologically obvious.
– uhoh
Jun 24, 2023 at 1:58
• @jlipinski thanks for the edit, it adds the relevant context for me Jun 26, 2023 at 10:30
• "I don't want to use straight mathematical approach" - please clarify - it is clearly a mathematical operation. Are you looking for the most efficient/streamlined/accurate approach (in which case it is a computational issue and not on topic). Or, are you just looking for a pre-written python function that does it (in which case it probably is on topic). Can you also define $\delta_0$? Jun 26, 2023 at 11:55
• @ProfRob I don't want to just paste this formula into wolfram and inverse it, resulting inversed function would be quite a mess. I'm looking for a function or commonly used formula for obtaining pixel coordinates knowing distortion coefficents $a_i$. $\delta_0$ is declination of corresponding center pixel, $\Delta\alpha$ is $\alpha - \alpha_0$ where $\alpha_0$ is rectascense of center pixel. Jun 26, 2023 at 12:36
• @jlipinski, so do you have a FITS header with the SIP coefficients? Jun 27, 2023 at 17:40

There is a fundamental mathematical problem that on paper is devastating and in practice irrelevant: there might not be any inverse, since clearly the $$\eta=f(x,y)$$ polynomial equation can have multiple solutions depending on the coefficients. Now, I bet the actual values of the coefficients merely generate a nonlinear but monotonous function that has an unique inverse.

The way I would solve this is to approximate that inverse function using linear interpolation at equidistant grid points, and then just get an interpolated answer. If the function is well-behaved over the grid the approximation error (estimable using the mean value theorem) will be small, and you can select the grid size so that error is smaller than the error due to the polynomial model.

This is how I would do it in Matlab: Let X be an array of the x values and Y the array of the y values, calculate eta=f(X,Y) as an array of eta values and delta similarly, and then just use the 2D linear interpolation x = interp2(eta, delta, X, etaq, deltaq) to find the x value corresponding to the query point etaq, deltaq.

• Yes, this is one possibility of approaching this problem, unfortunately it can be quite computational heavy. For 4k images there are over 8 milion pixels, that's a lot of etas and xis to calculate, but I agree, this method works. I'd prefer to wait few days for maybe some more answers. I know that astrometry engines (like astrometry.net) during plate solving are calculating coefficents for transforming both ways, unfortunately I have only this plate model, without access to original image. Jun 26, 2023 at 13:34
• Instead of interpolating (which generates a brand new (but small) polynomial each time you call it, one could also use the grid to fit a new inverse bi-quartic polynomial. It may be faster (or not) depending on floating point speed vs cache size/page fault speed.
– uhoh
Jun 26, 2023 at 18:22

Thank you all for your solutions and suggestions. During testing of Anders Sandberg solutions I've thought of another idea. I can calculate $$\xi$$ and $$\eta$$ directly from equatorial coordinates and thent I approximate X,Y using only linear coefficents of model (as there are easy to inverse), so:

$$\xi = a_1 + a_2x + a_3y \\ \eta = b_1 + b_2x + b_3y \\ x = \frac{\xi - (a_1 +a_3/b_3 *(\eta -b1))}{a_2-a_3 / b_3 * b_2} \\ y = \frac{\eta - b_1 - b_2 * x}{b_3}$$

Next I'm converting newly acquired X,Y again to equatorial using my model with every coefficent and comparing with queried RA, DEC equatorial coordinates. If the diffrence is small enough I accept those X,Y. If not, I'm using those new RA, DEC coordinates to approximate new X, Y, I'm calculating $$\Delta X = X_{new} - X_{old}$$ and $$X = X_{old} - \Delta X$$. Again I'm converting X,Y to RA DEC using my model and testing if the difference is small enough. It usually takes around one-two iterations to converge within 0.001 arcsec, but I think number of iterations may differ with more distorted images.

Anders Sandberg answer worked for me as well, but this one is much faster, as I only need to convert few points, not whole image.

• Sounds like a good approach. This is basically the secant method iterated a few times, for a function that is nearly linear. Jun 28, 2023 at 10:09
• Congratulations! Even though there are several good answers, when the OP solves the problem themselves it's great news!
– uhoh
Jul 2, 2023 at 9:12

If you have your SIP coefficients defined in a FITS header, this should be as simple as using the world_to_array_index_values() method of astropy.wcs.WCS.

The following is just an example since I don't have a fits file with SIP coefficients to test at the moment

import astropy.io.fits
import astropy.wcs

hdu_list = astropy.io.fits.open("myfile.fits")
hdu = hdu_list[0]