Creating a Milky Way Panorama on the Computer

Axel Mellinger
axm@rz.uni-potsdam.de

Feb 4, 2001

1  Introduction

The alternative approach is to take a series of pictures over the course of a year (possibly from different locations in order to capture both the northern and southern parts of the Milky Way), and later assemble them to a big panorama.

2  Image distortion

Let us look at a specific example to analyze the problems we are facing when combining several wide-angle pictures into a mosaic. Figure 1 shows two images of the Milky Way in the region between Canis Major and Crux with a generous overlap. The enlargements of the rectangular regions around the open cluster NGC 3114 reveal that the image scale near the edge is about 10-15% larger than at the center. Obviously, with such a large variation in image scale, we cannot seamlessly merge the two images.

The non-uniform image scale is the result of the camera lens projecting the celestial sphere onto planar film, as shown in Fig. 2. The three stars are separated by equal angles in the sky, but their linear separation in the film plane is larger for the pair that is further away from the optical axis. If we denote the angle between the star and the optical axis by s and the corresponding distance by s, the mathematical relationship between the two quantities is

s = f tans ,

(1)

where f is the focal length. Only the center of the image is free from distortion, since tans » s for small s.

3  Image transformation

The celestial coordinates (L,B) are related to the pixel coordinates (x,y) in the film plane through a unique mapping function F, i. e.

(x,y) = F(L,B)

(2)

If we knew this mapping function, we could take an arbitrary point in a Cartesian (rectangular), galactic coordinate frame, calculate the corresponding pixel (x,y) in the distorted raw image, and thus create a distortion-free image point by point. This is shown in Fig. 3: in the raw image (left), lines of constant galactic longitude and latitude are distorted into curves, whereas the transformed image (right) can be seamlessly merged with adjacent images.

We already know the mapping function in one dimension (cf. Eq. (1)). Unfortunately, reality is more complex; to work out the equations for a 2-dimensional film plane we have to resort to spherical trigonometry.

Figure 4 shows a view of the sky through the film plane. The position of the star is given by (x_s, y_s) in the film plane and by its galactic longitude L_s and latitude B_s in the sky. Likewise, the optical center is at (x_c,y_c) and (L_c,B_c), respectively. If we connect both the star and the image center to the galactic pole along lines of constant galactic longitude, we obtain a triangle where each side is part of a so-called great circle, i. e. a circle with the observer at its center. From a textbook on spherical trigonometry we obtain the relations for the position angle a and the angular distance s between the star and the optical axis:

coss = sinBs sinBc + cosBs cosBc cos(Ls-Lc) (3) sinscosa = sinBs cosBc - cosBs sinBccos(Ls-Lc) (4)

When solving for s and a, we have to pay special attention to their sign. While s is always positive, a has the same sign as DL. Now that we know s, we can calculate the distance s (in mm or pixels) between the star and the image center, using Eq. (1). Finally, x_s and y_s are obtained from s and a via

xs = xc - s sin(a-b) (5) and    ys = yc - s cos(a-b) , (6)

where x_c and y_c are the pixel coordinates of the image center and b is the angle by which the camera is rotated with respect to the galactic plane (b is positive for counter-clockwise rotation).

Unfortunately we cannot start to plug numbers into Eqs. (3)-(6) right away, as they contain a number of parameters that we do not accurately know. First, the galactic coordinates (L_c, B_c) of the optical axis can be roughly estimated using a star atlas, but for our purposes we need more precise numbers. Next, we need the position (x_c, y_c) of the optical image center on the film. It is not safe to assume that it coincides with the physical image center, since the image may have been shifted in the scanning process. Finally, we need accurate numbers for the focal length f of the camera lens and the camera rotation angle b.

The standard technique to determine the six parameters (L_c, B_c, x_c, y_c, f, b) is by means of a least squares fit. It is beyond the scope of this article to give a detailed introduction into least squares fitting procedures, but the basic idea is easily explained. As a first step, we pick n reference stars, where n is equal or greater than the number of parameters we try to determine, and measure their film and celestial coordinates (x_i, y_i,L_i, B_i), i = 1¼n. We also need a reasonably good initial guess for each parameter. With these initial values, we transform the (L_i, B_i) coordinates of our reference stars into film coordinates ([^(x)]_i, [^(y)]_i). The distance s_i between the computed and actual position of a reference star is given by

si2 = (xi- ^ x   i  )2 + (yi- ^ y   i  )2

(7)

The six parameters are now varied to achieve the ``best'' fit, i. e. the one that minimizes the sum over all s_i². A good introduction into least squares fitting techniques is given in Ref. [2].

Are we finally done? Not quite yet - our algorithm assumes a perfect lens, free of pincushion or barrel distortion, whereas most wide-angle lenses show at least moderate distortions. Other factors that may further distort the image are film shrinkage during processing and atmospheric refraction. A test calculation using an image taken with a 28 mm wide-angle lens and digitized to a size of 1500×1000 pixels showed an average error between the computed and actual film coordinates of 4 pixels, which is big enough to prevent accurate registration of the images. The quick remedy is to ``warp'' (or rather, de-warp) the image with a third-order polynomial transform, i. e.

x¢ = a00 + a10 x + a01 y + a20 x2 + a11 xy + a02 y2 + a30 x3 + a21 x2y + a12 xy2 + a03 y3 (8) y¢ = b00 + b10 x + b01 y + b20 x2 + b11 xy + b02 y2 + b30 x3 + b21 x2y + b12 xy2 + b03 y3  . (9)

The 20 coefficients a_ij and b_ij are again determined via a least squares fit. At least 10 reference stars are needed to uniquely specify the coefficients, but selecting 15 to 20 stars usually improves the stability of the calculation. After this ``warping'' step, the residual error was reduced to less than 1 pixel.

The complete algorithm has been coded into a PASCAL program which takes as input file a list containing the pixel coordinates, right ascension and declination of the reference stars as well as some parameters specifying the size of the transformed image. Conversion between equatorial (right ascension/declination) and galactic coordinates is done automatically.

4  Processing steps

1. Removal of dust specs and scratches
2. Reduction of background graininess.
   See http://canopus.physik.uni-potsdam.de/~axm/bgsmooth.html for details.
3. Flat-fielding. For a single image, a certain amount of vignetting is usually tolerable, but when several images are to be stitched together, special care must be taken to ensure a uniform background brightness. A number of anti-vignetting techniques have been described, e. g.
   • http://www.astropix.com/HTML/J_DIGIT/VIGNET.HTM
   • http://www.nanospace.com/~gottsch/vignet_en.htm
   • http://www.jps.net/jrfcomet/tech/vin.htm

4. Adjustment of color balance, background brightness and contrast. Only a rough adjustment is needed at this stage; final retouching can be done when combining the images in step 7.
5. Selection of reference stars. Aside from an image processing program that allows zooming in on individual stars, we need either a good star atlas and catalog or a planetarium program capable of displaying stars down to the limiting magnitude of the photograph. Wide-angle photographs typically go down to 9^th or 10^th magnitude. For Unix operating systems, XEphem is an excellent program which has support for the PPM catalog (limiting magnitude approx. 9^m) and the Hubble Guide Star Catalog (limiting magnitude approx. 14^m)
6. Image transformation to a galactic coordinate frame, as described in section 3.
7. Combining the transformed images. Most of the popular image processing programs, such as Adobe Photoshop or GIMP provide facilities for cutting and pasting images. To avoid visible edges, the selection should be feathered with a margin of at least 100 pixels.

Example

5  The final image

As described in section 3, the native map projection of the final image is a cylindrical equidistant projection. While allowing easy stitching of the individual frames, this projection has the disadvantage of showing substantial distortions in the form of elongated stars near the galactic poles. Therefore, two new versions of the image were created: A galacticAitoff projection and an azimuthal equidistant projection with two separate images for the northern and southern hemisphere.

The All-Sky panorama images can be seen at

http://home.arcor-online.de/axel.mellinger/allsky.html

A  Exposure Table

Sites

1. White Mountain Research Station (Barcroft Faciliy), California, USA (trip made possible by Tri-Valley Stargazers)
2. White Mountain Range (Grandview Campground), California, USA
3. Cederberg Observatory, Western Cape Province, South Africa
4. South African Astronomical Observatory, Sutherland, Northern Cape Province, South Africa
5. Swartberg Pass, Western Cape Province, South Africa
6. Organ Pipe Cactus National Monument, Arizona, USA
7. near Wünsdorf, Brandenburg, Germany
8. near Altes Lager, Brandenburg, Germany

References

[1]

J. Meeus, Astronomical Algorithms, Willmann-Bell, Inc.

[2]

W. H. Press, B. P.Flannery, S. A. Teukolsky and W. T. Vetterling, Numerical Recipes in PASCAL, Cambridge University Press, 1990

[3]

D. di Cicco, There's No Place Like Home, Sky and Telescope, Nov. 1999, pp. 137-140

[4]

A. Mellinger, There's No Place Like Home, Life Magazine, 2000 Alfred Eisenstaedt Awards for Magazine Photography, Spring 2000, pp. 34-35

[5]

A. Mellinger, Die Milchstraße im Computer, Sterne und Weltraum 39, Feb./March 2000, pp. 174-179

[6]

A. Mellinger, Un Panorama Galattico, Coelum 23, Oct. 1999, pp. 55-62