Pixel-to-RA/DEC mapping in digitized astrophotography

Question

I have a 1443x998 picture of the stars (taken w/ a 35mm camera and then scanned) with the following stars at the following pixel locations:

Altair x=782, y=532 [19h50m46.9990s RA, +08 52'05.959'' DEC] 
Sualocin, x=311, y=146 [20h 39m 38.287s +15 54'43.49'' DEC] 
Denebokab, x=1023, y=815 [19h25m29.9005s +03 06' 53.191'' DEC] 

What mathematical function converts pixel location to RA/DEC and vice versa? Notes:

• Bright stars are blobs in the picture; the coordinates above are roughly the center of the blob, but may be off by +-2 pixels.

• I know I can rotate the celestial sphere so that the center of my picture has polar coordinates 0,0. So the real question is "how to find this rotation" (but see next point).

• If elevation/azimuth were linear in pictures, this would be easy(er), but they're not: Measuring angular distance with photographs

• I can provide pixel locations of more stars if that helps. I believe 3 should be sufficient, but I could be wrong.

• I tried to choose 3 stars that were "spread out" across the picture (because I think that reduces error, not sure), but am not sure I succeeded.

• I'm doing this for several pictures and would like a general method.

• Doing this will help me identify fainter stars/Messier-objects/etc in the picture.

• I'm sure lots of astrophotographers want to do this, but haven't found any existing software that does this.

EDIT: Thanks, whuber! The gnomonic projection is what I was missing. I'd already done this assuming a linear transformation:

(* convert RA/DEC to xyz coords on celestial psuedo-sphere of radius 1 *) 
radecxyz[ra_,dec_] = 
{Cos[ra/12*Pi]*Cos[dec/180*Pi],Sin[ra/12*Pi]*Cos[dec/180*Pi],Sin[dec/180*Pi]}; 

(* I no longer have any idea how this works *) 
astrosolve[x_,y_,z_,xwid_,ywid_] := Module[{a,m,ans,nullans}, 
m=Array[a,{2,3}]; 
temp=Solve[{ 
m.radecxyz[x[[1]],x[[2]]]=={x[[3]]-xwid/2,x[[4]]-ywid/2}, 
m.radecxyz[y[[1]],y[[2]]]=={y[[3]]-xwid/2,y[[4]]-ywid/2}, 
m.radecxyz[z[[1]],z[[2]]]=={z[[3]]-xwid/2,z[[4]]-ywid/2} 
}]; 
ans = m /. Flatten[temp]; 
nullans=Flatten[NullSpace[ans]]; 
If[nullans.radecxyz[x[[1]],x[[2]]]<0,nullans=-nullans]; 
Return[{ans,nullans}]; 
]; 

where x, y, and z were each 4-element lists consisting of a stars RA, declination, x-coordinate on picture, and y-coordinate on picture. xwid and ywid are the width and height of the picture. In this case:

astrosolve[ 
 {19.8463886110, 8.8683219443, 782, 532}, 
 {20.6606352777, 15.9120805555, 311, 146}, 
 {19.4249723610, 3.1147752777, 1023, 815}, 
 1443, 998] 

{ 
 {{-2250.51, -1182.52, 385.689},  {-166.12, -543.746, -2376.73}},  
 {0.480698, -0.861509, 0.163497} 
} 

Now, referring to "{-2250.51, -1182.52, 385.689}" as $frow, "{-166.12, -543.746, -2376.73}" as $srow, and "{0.480698, -0.861509, 0.163497}" as $null, this PHP subroutine translates RA/DEC to xy coordinates:

# radecxy(ra,dec): converts ra/dec to x,y using a quasi-linear transformation 

function radecxy($ra,$dec) { 
    global $null,$frow,$srow,$xwid,$ywid; 
    list($x,$y,$z)=array(cos($dec)*cos($ra),cos($dec)*sin($ra),sin($dec)); 

    $dotprod=$null[0]*$x+$null[1]*$y+$null[2]*$z; 
    if ($dotprod<0) {return(array(-1,-1));}

 list($fx,$fy)  = array($frow[0]*$x+$frow[1]*$y+$frow[2]*$z,$srow[0]*$x+$srow[1]*$y+$srow[2]*$z); 
    $fx+=$xwid/2; 
    $fy+=$ywid/2; 
    if ($fx<0 || $fy<0 || $fx>$xwid || $fy>$ywid) { 
        return(array(-1,-1)); 
    } else { 
        return(array($fx,$fy)); 
    } 
} 

Sadly, I no longer have any idea why this works, but using it + adding known star positions yields tolerable results (use "view image" to see it full size):

However, as you can see, the results aren't perfect, convincing me that a linear transformation wasn't the right answer. I think gnomonic might be the grail I was seeking.

Answer 1

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.

Answer 2

To do this with the same degree of precision that professional astronomers do it would indeed be difficult. It would require you to have extremely precise characterization of the distortions produced by your lens, and the imperfections in your camera's sensor. However, you probably don't need that degree of accuracy. It should be sufficient for you to assume that your lens does not introduce large amounts of distortion (which is a good assumption for a quality lens) and that your cameras sensor is pretty close to a perfectly regular grid (which is a very good assumption for even a cheap camera).

All that remains is to work out the coordinate transformation that describes the orientation of the camera i.e. the direction it was pointed and the degree to which it was rotated.

What you are looking for then, is called an affine transformation, or an affine map. Which is just a fancy name for a matrix by which you would multiply your pixel coordinates to obtain your astronomical coordinates. In the case of an affine map, this transformation can include any degree of rotation, scale, shear, and translation.

The meaning of the rotation component is pretty obvious. The scale factor simply describes how much of the sky is covered by each pixel in terms of RA/Dec. Shear is a transformation that would make the image of a rectangle become a parallelogram, but there shouldn't be any of this effect in an image of objects at infinity (like stars). Lastly, the translation component simple adds an offset to account for the fact that the (x = 0, y = 0) pixel in your image probably does not correspond to (RA = 0, Dec = 0).

Because you have 3 reference stars in your image, you have enough information to calculate the relationship between your pixel coordinates and the RA/Dec you are looking for. This would be done by linear least squares fit (not non-linear least squares as mentioned above) to determine the values of the matrix components that best match your pixel coordinates to the known RA/Dec of the reference stars. Once the matrix is established you can then apply it to the pixel coordinates of other stars to get their RA/Dec.

While I could do this relatively easily, I'm unfortunately unsure how to help you do it. It would involve some mathematical skill that is a bit beyond the scope of photo.SE. I'm an optical engineer, but I'm not much of a photographer; the software I would use for this is designed for engineers to do heavy-duty numerical computation, and isn't really a photographic tool at all. There may be ways to do this using software packages geared towards photographers, but I don't know about them.