f



Least square fitting

Hi all,

I have a problem with a rather complicated function depending on four
parameters which I try to find using least square fitting and I don't
know exactly how to do it.

The basic problem is the following:
I have an astronomical image of a star field and try to relate the sky
coordinates (right ascension, declination) of the stars to the pixel
coordinates (x, y).

The function to relate the two depends on the not accurately known
parameters focal length of the lens (f), the rotation of the field of
view with respect of the north direction (beta) and the center sky
coordinates of the image (bc, lc).
I know the sky coordinates of certain stars accurately (bs, ls) as
well as their pixel coordinates (xs, ys) and the center pixel
coordinates of the image (xc, yc).

Generally, the function is as follows (equivalent for ys', but the
functions f1-f5 are slightly different):

xs' = f1(xc, f2(f, f3(ls, lc, bs, bc)), f4(ls, lc, bs, bc, f5(ls, lc,
bs, bc)), beta)

The functions f1-f5 look rather complicated and contain many sines,
cosines and acos etc.

I then did the following as I didn't know better:
First I minimized all this with respect of beta, then lc and bc and
after that f starting with some initial values for the four
parameters. Then I started all over again, beta, lc, bc, f (in a loop)
until the difference between the known coordinates (xs, ys) and the
calculated ones (xs', ys') reached a minimum.

However, the minimum deviation reached at the end of the loop pretty
much depends on the initially chosen values of beta, lc, bc and f.

There certainly is a more optimal way to do the least square fitting
than the one that I have chosen (partial derivatives etc.), but that
f1 depends on f2 and f4 and these again depend on f3 and f5 is giving
me some unsolvable problems.

How should I proceed? What would be a good method to minimize this?

Thanks, Michael


The main problem now is for me, that f1 depends on f2 and f4 and that
f2 and f4 themselves depend on f3 and f5.

0
5/20/2008 9:41:32 AM
comp.lang.idl-pvwave 12260 articles. 4 followers. Post Follow

7 Replies
448 Views

Similar Articles

[PageSpeed] 5

On May 19, 11:41 pm, MichaelT <michael.theus...@googlemail.com> wrote:
> Hi all,
>
> I have a problem with a rather complicated function depending on four
> parameters which I try to find using least square fitting and I don't
> know exactly how to do it.
>
> The basic problem is the following:
> I have an astronomical image of a star field and try to relate the sky
> coordinates (right ascension, declination) of the stars to the pixel
> coordinates (x, y).
>
> The function to relate the two depends on the not accurately known
> parameters focal length of the lens (f), the rotation of the field of
> view with respect of the north direction (beta) and the center sky
> coordinates of the image (bc, lc).
> I know the sky coordinates of certain stars accurately (bs, ls) as
> well as their pixel coordinates (xs, ys) and the center pixel
> coordinates of the image (xc, yc).
>
> Generally, the function is as follows (equivalent for ys', but the
> functions f1-f5 are slightly different):
>
> xs' = f1(xc, f2(f, f3(ls, lc, bs, bc)), f4(ls, lc, bs, bc, f5(ls, lc,
> bs, bc)), beta)
>
> The functions f1-f5 look rather complicated and contain many sines,
> cosines and acos etc.
>
> I then did the following as I didn't know better:
> First I minimized all this with respect of beta, then lc and bc and
> after that f starting with some initial values for the four
> parameters. Then I started all over again, beta, lc, bc, f (in a loop)
> until the difference between the known coordinates (xs, ys) and the
> calculated ones (xs', ys') reached a minimum.
>
> However, the minimum deviation reached at the end of the loop pretty
> much depends on the initially chosen values of beta, lc, bc and f.
>
> There certainly is a more optimal way to do the least square fitting
> than the one that I have chosen (partial derivatives etc.), but that
> f1 depends on f2 and f4 and these again depend on f3 and f5 is giving
> me some unsolvable problems.
>
> How should I proceed? What would be a good method to minimize this?
>
> Thanks, Michael
>
> The main problem now is for me, that f1 depends on f2 and f4 and that
> f2 and f4 themselves depend on f3 and f5.

How big is the field of view? If it is small enough (and the optics
consistent enough) that you don't need to worry about the pixel scale
(number of arcsecs/pixel) changing over the image, than the problem
seems easier, since most of the parameters to solve for ought to be
independent. I don't know how statistically kosher this is, but here
are some thoughts:

Assuming that the pixel scale p in arcsec/pixel is a constant, then
the RA,DEC (a,d) of stars are related to their pixel coords (x,y) by

a=ac+((x-xc)*cos(beta)-(y-yc)*sin(beta))*p
d=dc+((x-xc)*sin(beta)-(y-yc)*cos(beta))*p

where ac,dc,xc,yc are the sky and pixel coordinates of the center
pixel. Beta is the angle of rotation of the image y axis with respect
to north. I guess also here I've assumed that RA increases to the
right (backwards from convention) and that the image is near the
celestial equator so that 1 second of right ascension = 1 second of
declination. Both of these are easy to generalize (switch some of the
signs for East pointing left, and add in a cos(dec) term for off-
equator correction). Neither of these should affect the following
analysis.

For every (a,d) (x,y) pair, compute the pixel distance and angular
separation. A simple linear regression then determines p (an easier
but sloppier method would be to calculate p for each pair and take a
weighted average). p is independent of all of the other parameters, so
(to me) it seems like you could do this seperately.

Now that you know p, for a bunch of pairs of stars, compute the x
displacement in pixels and the RA displacement in arcsec. Convert the
x displacement in pixels to an angular distance using p. Now use the
equations above to solve for beta. (you can also use the y coordinates
and/or decs just as easily).

Finally, determining ac and dc is one more linear regression.

Is this the sort of thing you are trying to do?

Chris

0
cnb4ster (16)
5/20/2008 11:16:37 AM
Thanks for your thoughts, Chris!

But I think the field of view is not small enough to do it this way
(about 3.5=B0 x 2.5=B0 and larger). The focal length is shorter than 400
mm. Otherwise, it could very well work the way you described it, I
think.

These are the things I'd like to do in the end:
I have several images of the same object, but from different nights
taken with an amateur instrument. So the fov is slightly different
each time as well as the rotation angle (it's not a stationary scope).
Knowing all the parameters would enable me to exactly overlay the
images to increase the signal-to-noise ratio. I know that there are
various free programs (non-IDL) that can already do that, but they all
have some shortcomings. So I tried to do something myself.
Also, you could then use some neat routines from the IDL astronomy
library (http://idlastro.gsfc.nasa.gov/) and automatically generate an
overlay from the USNO catalog to determine the limiting magnitude or
to label objects (e.g. queryvizier and querysimbad).
There are many other things you could do if you knew the pixel-to-
coordinate conversion.

Michael
0
5/20/2008 2:06:48 PM
On May 20, 5:41 am, MichaelT <michael.theus...@googlemail.com> wrote:

> The basic problem is the following:
> I have an astronomical image of a star field and try to relate the sky
> coordinates (right ascension, declination) of the stars to the pixel
> coordinates (x, y).
>
I have performed a lot of this type of astrometry, and your problem
seems very strange to me.

> The function to relate the two depends on the not accurately known
> parameters focal length of the lens (f), the rotation of the field of
> view with respect of the north direction (beta) and the center sky
> coordinates of the image (bc, lc).

A minor note is that if you really have such a large distortions then
a single beta does not make sense -- the rotation from North will vary
slightly across the image.


>
> Generally, the function is as follows (equivalent for ys', but the
> functions f1-f5 are slightly different):
>
> xs' = f1(xc, f2(f, f3(ls, lc, bs, bc)), f4(ls, lc, bs, bc, f5(ls, lc,
> bs, bc)), beta)
>
> The functions f1-f5 look rather complicated and contain many sines,
> cosines and acos etc.
>
> I then did the following as I didn't know better:
> First I minimized all this with respect of beta, then lc and bc and
> after that f starting with some initial values for the four
> parameters. Then I started all over again, beta, lc, bc, f (in a loop)
> until the difference between the known coordinates (xs, ys) and the
> calculated ones (xs', ys') reached a minimum.
>

I have no idea what your functions f1, f2, f3, f4 ,f5 might be, but do
you really care about them?    All you want is a function that maps
the known coordinates to the calculated ones.    Often what one does
is to find a quadratic or cubic function that will linearize the x,y
coordinates (i.e. so that they line up with RA, Dec)

xp = x + axy + by^2 + c*x^2 + ..
yp = y + dxy + ey^2 + f*x^2

You then use least squares to determine the a,b,c.. coefficients.
This is the 'SIP' convention discussed in
http://ssc.spitzer.caltech.edu/postbcd/doc/shupeADASS.pdf and used
e.g. by the astrometry.net software for handling distortions.    Of
course, it would be nice to know the true functional form of the
distortions but whatever functions you are using don't seem very
useful.     You can think of the quadratic or cubic equations as a
Taylor series approximation to whatever the true functional form is.


--Wayne
0
wlandsman (223)
5/20/2008 5:28:32 PM
> I have performed a lot of this type of astrometry, and your problem
> seems very strange to me.
The way I proceed is described here:
http://home.arcor-online.de/axel.mellinger/mwpan_web/mwpan_web.html
It uses the true functional form of the conversion (ra, dec) -> (x,y).
So I am not so sure if we are both talking about the same thing here.
I am not after the distortions, yet (These seem to be small. The
calculated x',y' positions deviate by only 0.5 pixels, on average,
from the real ones - when the algorithm converges).

> you really care about them?    All you want is a function that maps
> the known coordinates to the calculated ones.    Often what one does

What I want is a function that converts (ra, dec) into pixel
coordinates (x, y). So I am not so sure if that is what you describe
in the following. What you describe looks more like (x', y') -> (x, y)
to me.

> is to find a quadratic or cubic function that will linearize the x,y
> coordinates (i.e. so that they line up with RA, Dec)
>
> xp = x + axy + by^2 + c*x^2 + ..
> yp = y + dxy + ey^2 + f*x^2
>
> You then use least squares to determine the a,b,c.. coefficients.
> This is the 'SIP' convention discussed inhttp://ssc.spitzer.caltech.edu/postbcd/doc/shupeADASS.pdfand used
> e.g. by the astrometry.net software for handling distortions.

The webpage I posted above also discusses the elimination of the
distortions and uses a cubic function of the form given by you above
and as discussed in the paper.

Can I also use the cubic/quadratic function to map (ra, dec) to (x, y)
or vice versa?

This is my first attempt as an amateur astronomer to deal with these
things :-)
Many thanks for your help.

Michael
0
5/20/2008 6:44:40 PM
MichaelT <michael.theusner@googlemail.com> wrote:
>  Hi all,
> 
>  I have a problem with a rather complicated function depending on four
>  parameters which I try to find using least square fitting and I don't
>  know exactly how to do it.
> 
>  The basic problem is the following:
>  I have an astronomical image of a star field and try to relate the sky
>  coordinates (right ascension, declination) of the stars to the pixel
>  coordinates (x, y).

So if you truly wish to know how to properly solve this particular
mathematical problem, I will leave that to other wiser experts to 
answer. 

If, however, what you really want is just to know the astrometric
solution for your data, you may be interested to know about
Astrometry.net, a research project by David Hogg at NYU and
collaborators which intends to solve this particular problem once and
for all, for everyone. You feed their web service arbitrary
astronomical images, and it hands back the astrometric solution. I've
not tried it for images as wide-field as you have, but I suspect it
will work just fine.  See http://astrometry.net/ for more info.

  - Marshall
0
mperrin (13)
5/20/2008 9:48:34 PM
Marshall Perrin wrote:
> MichaelT <michael.theusner@googlemail.com> wrote:
>>  Hi all,
>>
>>  I have a problem with a rather complicated function depending on four
>>  parameters which I try to find using least square fitting and I don't
>>  know exactly how to do it.
>>
>>  The basic problem is the following:
>>  I have an astronomical image of a star field and try to relate the sky
>>  coordinates (right ascension, declination) of the stars to the pixel
>>  coordinates (x, y).
> 
> So if you truly wish to know how to properly solve this particular
> mathematical problem, I will leave that to other wiser experts to 
> answer. 
> 
> If, however, what you really want is just to know the astrometric
> solution for your data, you may be interested to know about
> Astrometry.net, a research project by David Hogg at NYU and
> collaborators which intends to solve this particular problem once and
> for all, for everyone. You feed their web service arbitrary
> astronomical images, and it hands back the astrometric solution. I've
> not tried it for images as wide-field as you have, but I suspect it
> will work just fine.  See http://astrometry.net/ for more info.
> 
>   - Marshall
I agree with Marshall. Astrometry.net offers a good service and the
code works well if the FOV is not too small (i.e. there are enough
stars to match). However, it does not take (yet) advantage of the WCS
info provided by the FITS header. For your purpose the TERAPIX
procedures SCAMP and SWARP can be used. The former establishes
the astrometry (including field distortion) while the latter does the image
transformation, i.e. creates a non-distorted output mosaic. You might
want to have a look at the TERAPIX website for more information.

Regards,

	Bringfried
0
stecklum (65)
5/21/2008 7:47:53 AM
Thanks Marshall and Bringfried,

I will try it out. Actually, I really only need to know the conversion
parameters. Everything else can then be done using those.

It would certainly be nice to do it all yourself, but, selecting the
reference stars is always a hassle.

Michael
0
5/21/2008 8:45:56 AM
Reply: