<p>There is a fundamental mathematical problem that <em>on paper is devastating and in practice irrelevant</em>: there might not be any inverse, since clearly the <span class="math-container">$\eta=f(x,y)$</span> 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.</p>
<p>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.</p>
<p>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.</p>

CC BY-SA 4.0
2023-06-26T13:34:49.343Z
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.
jlipinski
CC BY-SA 4.0
2023-06-26T18:22:07.717Z
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
CC BY-SA 4.0
2023-06-26T13:15:11.900Z
2023-06-26T13:15:11.900Z
Anders Sandberg
<p>Thank you all for your solutions and suggestions. During testing of Anders Sandberg solutions I've thought of another idea. I can calculate <span class="math-container">$\xi$</span> and <span class="math-container">$\eta$</span> directly from equatorial coordinates and thent I approximate X,Y using only linear coefficents of model (as there are easy to inverse), so:</p>
<p><span class="math-container">$
\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}
$</span></p>
<p>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 <span class="math-container">$\Delta X = X_{new} - X_{old}$</span> and <span class="math-container">$X = X_{old} - \Delta X$</span>. 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.</p>
<p>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.</p>

CC BY-SA 4.0
2023-06-28T10:09:28.323Z
Sounds like a good approach. This is basically the secant method iterated a few times, for a function that is nearly linear.
Anders Sandberg
CC BY-SA 4.0
2023-07-02T09:12:50.390Z
Congratulations! Even though there are several good answers, when the OP solves the problem themselves it's great news!
uhoh
CC BY-SA 4.0
2023-06-28T08:42:49.710Z
2023-06-28T08:42:49.710Z
jlipinski
<p>If you have your SIP coefficients defined in a FITS header, this should be as simple as using the <code>world_to_array_index_values()</code> method of <code>astropy.wcs.WCS</code>.</p>
<p>The following is just an example since I don't have a fits file with SIP coefficients to test at the moment</p>
<pre><code>import astropy.io.fits
import astropy.wcs

hdu_list = astropy.io.fits.open(&quot;myfile.fits&quot;)
hdu = hdu_list[0]
wcs = astropy.wcs.WCS(hdu.header)
x, y = wcs.world_to_array_index_values(xi, eta)
</code></pre>

CC BY-SA 4.0
2023-06-28T01:12:22.617Z
2023-06-28T01:12:22.617Z
Roy Smart
<p>I have a polynomial model that can transform pixel X,Y coordinates to standard <span class="math-container">$\xi, \eta$</span>. This model looks like this:</p>
<p><span class="math-container">$\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$</span></p>
<p>where I know every <span class="math-container">$a_i$</span> cooefficent for both <span class="math-container">$\xi$</span> and <span class="math-container">$\eta$</span>. I also know equatorial coordinates of center pixel so I can use those formulas to get equatorial coordinates of given X, Y pixel:</p>
<p><span class="math-container">$
\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}}
$</span></p>
<p>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?</p>
<p>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 <span class="math-container">$a_i$</span> 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).</p>
<p>EDIT2:
Sorry for misunderstanding. I don't have corresponding image or valid FITS header. Only thing I have got is SIP (the <span class="math-container">$a_i$</span> coefficetns and <span class="math-container">$\alpha_0$</span> and <span class="math-container">$\delta_0$</span>) model for converting from pixel to standard.</p>

CC BY-SA 4.0
2023-06-24T01:58:20.910Z
@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
CC BY-SA 4.0
2023-06-26T10:30:11.760Z
@jlipinski thanks for the edit, it adds the relevant context for me
planetmaker
CC BY-SA 4.0
2023-06-26T11:55:08.943Z
"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$?
ProfRob
CC BY-SA 4.0
2023-06-26T12:36:32.200Z
@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.
jlipinski
CC BY-SA 4.0
2023-06-27T17:40:48.953Z
@jlipinski, so do you have a FITS header with the SIP coefficients?
Roy Smart
CC BY-SA 4.0
2023-06-28T07:15:25.933Z
Please add that you don't have the relevant image. It's kind-of important.
ProfRob
<coordinate><astrometry><frame-of-reference><wcs><coordinate-systems>
Reverse polynomial model for reverse astrometry
jlipinski
# Reverse polynomial model for reverse astrometry

(Asked by: username_0 on 2023-06-23)

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.

**Tags:** coordinate, astrometry, frame-of-reference, wcs, coordinate-systems
> **Comments:**
> **username_1**: @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.  
> **username_2**: @jlipinski thanks for the edit, it adds the relevant context for me  
> **username_3**: "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$?  
> **username_0**: @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.  
> **username_4**: @jlipinski, so do you have a FITS header with the SIP coefficients?  
> **username_3**: Please add that you don't have the relevant image. It's kind-of important.  

---

## 3 answers

### Answer by username_5 on 2023-06-26. Score: 3 | Accepted answer

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.

> **Comments:**
> **username_0**: 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.  
> **username_1**: 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.  

---

### Answer by username_0 on 2023-06-28. Score: 3

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.

> **Comments:**
> **username_5**: Sounds like a good approach. This is basically the secant method iterated a few times, for a function that is nearly linear.  
> **username_1**: Congratulations! Even though there are several good answers, when the OP solves the problem themselves it's great news!  

---

### Answer by username_4 on 2023-06-28. Score: 2

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]
wcs = astropy.wcs.WCS(hdu.header)
x, y = wcs.world_to_array_index_values(xi, eta)


---

