I added experimental support for astrometry to Munipack. It is primary designed for astrometry on CCD. So only gnomonical projection and affine transformation (shift and rotation) are supported yet. Both limitation can be freed in future releases.
The astrometry code is available in GoogleCode now.
The algorithm
I coded my own implementation of the widely-known "Astrography plate
measurement" (Smart & Green, Tagg).
On input:
The table with Right Ascension,Declination (from any catalogue) and corresponding measured rectangular x and y coordinates for every object.
The algorithm:
- Get center of a picture in pixels as half of both sizes.
- Estimate a center of projection in spherical coordinates by mean.
- Project spherical coordinates to rectangular ones.
- Estimate scale as mean ratio between distances of projected and measured distances of objects.
- Estimate angle of rotation between projected and measured coordinates by using of property of scalar product of vectors.
- Compute the transformation by some minimization method with starting values given by items 3,4 and 5.
- Check offsets between center of projected and measured coordinates. If one is too great, does better estimate of the center of projection and repeat from item 2, else finish.
Precise coordinates of center of projection, scale and position angle of image. Optionally, statistical parameters.
The described algorithm is relative general, so it can be easy modifies for another projections or more complex transformations of rectangular coordinates like distortions, pincushion, etc.
The algorithm is iterative due to fact that we need known the center of projection before we have calibrated image. Usually, only tree iterations are required.
The Robust Algorithm
The real-world usable algorithm needs to be robust. Small deviations can cause only small variances in output parameters. Outliners has only small fluency on the solution. The simple example of the use of robust method can be found in Launer & Wilkinson: Robustness in statistics: proceedings of a workshop (1979) or in Numerical Recipes. Munipack implements robust mean estimator in lib/statistics.f90 as rmean. Another use of the robust algorithm is on straight light fitting.
I applied the basic schema of robust algorithms to the above algorithm:
- Estimate median of absolute deviations (MAD) using of minimization of absolute deviation (Nelder-Mead).
- Use method of maximal likelihood to estimate the proper transformation.
Astrometry implements module in astrometry/astrometry.f90. The estimation of MAD is loop around nelmin. The maximum likelihood is loop around Minpack's hybrd.
The fitting is decomposed on tree single steps (estimate of center of projection, scale and affine transformation) so we works in parameter sub-spaces (we are not fitting all the parameters together). Numerical experiments gives me more robust behavior and more precise solutions than simultaneous fitting of all parameters. By my opinion, it is results of (non-)ignoring of their cross-correlations. It is interesting that QR decomposition does not solve the problem.
The algorithm is optimized on precision not on speed. So for a lot of objects may be slow.
The Low Precision Algorithm
The robust algorithm can be used only for many object (ideally for tens of stars and more). The absolute minimum is five stars. I'm supposing that sometimes will be required astrometry for less stars. In this case, I'm switching to a simple algorithm which uses median only to estimate of required parameters.
For two stars, the transformation can be determined, but is is not possible to estimate uncertainties of parameters. It is not possible to estimate scale or rotation for only single star.
Input/Output "protocols"
I choose an unusual way to set an input and get an output data. The stream input is a text file consists a set of lines (records) with the same structure as FITS headers. By the user point of view, the routine performs as a transformation filter. The text file on input is modified by replacing and adding new lines to the output text file.
The output structure is perfectly suitable to be directly writable to a FITS header in standard WCS format. We need an external utility to merge of the output lines to an existing FITS file. Munipack provides munifits tool for the task.
It is possible to add another type of output, but I think, there is no other widely used standard astrometric format. If the way will successful, I'll reimplement it for remaining utilities of Munipack.
The style can be little bit strange but absolutely the same as standard e-mail handling. An user interface (Thunderbird, elm, ..) creates well formatted input file and pass it to a SMTP daemon. The daemon (or its successors) delivers its (as copying filter) to a recipient and some client decode its. Note that many Unix utilities uses the same idioms (sed, awk,...). The main advantage of the style is a simple modification of format, also back and forward compatibility. [http://catb.org/~esr/writings/taoup/].
The protocol is in detail described in included manual page and test example.
An example of calibrated image with UCAC2 displayed stars. Thx. Mr. Pavel. Large size.
No comments:
Post a Comment