I will outline a rigorous approach and indicate what software can help with it. Most of this will be tangential to the interests of the photography site, but because there are some useful insights that apply to any circumstance in which locations will be estimated from measurements on an image, this site seems a reasonable place for such an analysis.

Taking an image (with a lens that has been corrected for distortion) projects the celestial sphere through the focal point of the lens onto the plane of the sensor. This is an oblique aspect of a gnomonic projection.

Mathematically, the conversion from (RA, DEC) proceeds through a series of steps:

1. Convert (RA, DEC) to spherical coordinates. RA has to be converted from hours-minutes-seconds into degrees (or radians) and DEC has to be converted from degrees-minutes seconds into degrees (or radians), remembering that it is elevation above the plane, not angle from the north pole (which is the usual spherical coordinate convention). Both conversions are simple arithmetic.

2. Compute (x,y,z) coordinates for the spherical coordinates of the stars. This is a standard coordinate conversion (involving simple trigonometry).

3. Rotate the celestial sphere to align its poles with the lens axis. This is a linear transformation.

4. Rotate the celestial sphere around its poles to conform to the camera's orientation (another linear transformation).

5. Placing the imaging plane at a constant height z above the focal point, draw light rays from the stars at (x,y,z) through the focal point until they intercept the plane. (This is the gnomonic projection and, by its nature, it is projective and not linear.)

[In the figure, which is intended to be a planar cross-section through the axis of the lens,

• A is the focal point.
• Semicircle BCD is the visible part of the celestial sphere.
• AC points along the axis of the lens.
• E, F, and G are star locations.
• EE, FF, and GG are their corresponding locations on the (invisible) celestial sphere.
• E', F', and G' are their images on the sensor KL (so that EE', FF', and GG' are paths of light rays from the stars to the sensor).
• AD is the horizon from which declination is measured.
• Alpha is the declination of star E (or, equivalently, an angular coordinate of EE). Stars F and G have similar declinations (not shown).

Our task is to find the mathematical relation between the angular coordinates for E, F, and G--which are assumed known to high accuracy--such as alpha, and coordinates of their images E', F', and G', measured in pixels along the sensor. Once found, this relation can be inverted as described below to estimate angular coordinates of celestial objects from positions of their images on the sensor. Not shown, for simplicity, is the magnification of the lens. With a distortion-free lens this will have the effect of uniformly rescaling the coordinates of E', F', and G' relative to the center of the sensor.]

This procedure describes how the light gets from a star onto the sensor for a perfect simple lens. It involves these (unknown) parameters, which will need to be determined:

• Three angles in (3) and (4) describing the lens and camera orientation.

• A scale factor in (5) describing the combined effects of the sensor's size, distance from the focal point, and magnification of the lens.

Due to the projection (5), this is a complex, nonlinear transformation in general, but it has a definite mathematical description. If we let x = (RA, DEC) designate a star's position, let theta represent the four parameters for the imaging process, and let y = (column, row) represent the pixel coordinates, then we can abstractly but more simply write

y = f(x, theta).

Next--and this is highly important--we need to account for errors. The imaged stars are not in precise locations. Thus we have to include an error term in our formula and it's conventional (since about 1800) to model this error probabilistically. The new formula is

y = f(x, theta) + e

When the lens is distortion free, the expected value of e is 0 and its standard deviation (sigma) measures the typical size of error. It's reasonable to assume the e's are approximately normally distributed, with approximately equal standard deviations (which is not true, but for an initial analysis it's a reasonable assumption) and we can hope these errors are statistically independent of each other (which again is not true but it's a good starting assumption). This justifies a least squares solution using maximum likelihood. Up to a universal constant, whose value we don't need to know, the log probability of any particular observation (x,y) equals

-|f(x, theta) - y|^2 / (2 sigma^2) - 2 log(sigma).

(The absolute value bars denote Euclidean distance in the imaging plane, computed as usual with the Pythagorean Theorem.)

By virtue of the assumed independence of errors, the log probability of the set of data for an image is the sum of these values. This is the "log likelihood." The maximum likelihood (ML) estimates of the parameters theta and sigma (five numbers in all) are those values that maximize the log likelihood.

We can, and should, go further. The theory of ML also shows how to obtain confidence intervals for the estimates. Intuitively, the errors in our observations create a little uncertainty in the joint values of the angles, the scale factor, and the standard deviation. We need these values to estimate RA and DEC for any pixels in our image. By using uncertain values, which is unavoidable, we will get uncertain results. Furthermore, if we identify a pixel in our image by looking at a diffuse blob of light (scattered over approximately pi*sigma^2 pixels altogether), there will be additional uncertainty in the pixel coordinates. Collectively these two forms of uncertainty combine. This implies the net uncertainty in estimating the RA and DEC of any blob of light on the image is larger than you would guess.

Finally, when you're taking a measurement off the image and using that to estimate the true coordinates of a star or celestial object, you are doing inverse regression, which is a form of instrument calibration. Inverse regression is a procedure to account for the uncertainties I just described. Its output of course includes the estimated star coordinates for any blob of pixels on the image. It also includes a ring of coordinates around that estimate that are also consistent with the location of that blob. (This is a joint "inverse prediction interval" or a set of "fiducial limits" for the blob's RA and DEC.) In practice, if you consult a catalog of celestial objects, you can use this ring to search for all known objects that are consistent with the information in your image. Clearly this can be more valuable than a simplistic procedure that estimates--sometimes incorrectly--only one single set of coordinates.

In summary, what is needed here is software to

• Perform the nonlinear optimizing required by ML.

• Estimate standard errors in the estimates.

• Perform inverse regression.

Expertise with appropriate software, such Stata's ML command or Mathematica, is essential if you're coding this yourself.

Regardless of your expertise, here are some conclusions you can use in your imaging strategies:

• The accuracy of the image for getting a fix on any object can never be greater than the inherent imprecision in the imaging (as measured by sigma, the typical size of a point of light on the image).

• You can get close to this level of accuracy by identifying many known stars, not just three. This reduces the uncertainty in the sky-to-image transformation almost to zero if you have enough known stars located in the image.

• It is correct that you want the reference stars to be spread across the image. It is also crucial that they not be lined up (which, unfortunately, is the case with the three locations given in the question). If you can afford to locate only three stars, get them in a nice triangle. When the stars do line up, the statistical analysis indicates there is a huge uncertainty about locations in directions perpendicular to the line. In this particular example the estimated error (sigma) is hundreds of pixels wide. Having one more star to make a good triangle ought to reduce this error to one or two pixels.

Some parting thoughts:

• It is possible to detect and even correct lens aberrations by carrying out a more extensive statistical analysis. The idea is to plot deviations between expected and actual locations of the stars on the image. This is akin to "warping" or "georeferencing" map data. As a quick and dirty solution, you can press GIS or image processing software (such as ENVI) into service to georeference (or astroreference) any image. Such software usually does not carry out ML estimates of projective transformations, but it can do high-order polynomial approximations. An order-2 or order-3 polynomial transformation might do a good enough job, depending on your application.

• It is possible to improve accuracy by combining multiple images of the same objects.