Creating a Milky Way Panorama on the Computer

Axel Mellinger
axm@rz.uni-potsdam.de

Feb 4, 2001

1  Introduction

Looking at the Milky Way from a dark (preferably southern) site can be one of the most memorable experiences. Many astrophotographers feel the wish to capture the panoramic view on film. However, the large area covered by the Milky Way cannot be easily photographed on a single frame. Table 1 shows the field of view of some popular camera lenses. Even with a 28 mm wide-angle lens, we can capture a mere 65 degrees of the Milky Way. Of course, there are super-wide-angle and fish eye lenses, but in order to achieve the very wide field of view, some compromises in the optical design have to be made, and they tend to be rather slow.

focal length field of view
135 mm 15×10
50 mm 40×27
28 mm 65×46
Fish eye 180
Table 1: Field of view of some popular camera lenses

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

img51.jpg img36.jpg
 
img51mag.jpg img36mag.jpg

Figure 1: Milky Way between Crux and Canis Major. Even though both images were taken with the same 28 mm lens, the image scale of the yellow rectangles is different, as can be seen in the enlargements.

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.

1dproject.gif
Figure 2: Projection of the celestial sphere onto the film plane (top view).

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

In order to correct the image scale distortions discussed in the previous section, we need to work out an algorithm relating the position of a star on the film (measured either in mm or pixels) to its celestial coordinates. The position of each star in the sky can be described by a pair of coordinates, such as right ascension and declination. For a Milky Way panorama, the galactic longitude L and latitude B are the most suitable coordinates. Conversion formulae between the two systems are given by most textbooks on spherical or computational astronomy (see, e. g., [1]).

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.

distort.gif undistort.gif

Figure 3: Lines of equal galactic longitude and latitude on the film plane. Left: raw image, right: ideal (equidistant cylindrical) projection. The camera is assumed to be centered on the galactic equator (B = 0).

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.

2dproject.gif
Figure 4: Projection of the celestial sphere onto the film plane as seen from a place behind the camera.

Figure 4 shows a view of the sky through the film plane. The position of the star is given by (xs, ys) in the film plane and by its galactic longitude Ls and latitude Bs in the sky. Likewise, the optical center is at (xc,yc) and (Lc,Bc), 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, xs and ys are obtained from s and a via
xs
=
xc - s sin(a-b)
(5)
and    ys
=
yc - s cos(a-b) ,
(6)
where xc and yc 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 (Lc, Bc) 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 (xc, yc) 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 (Lc, Bc, xc, yc, 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 (xi, yi,Li, Bi), i = 1n. We also need a reasonably good initial guess for each parameter. With these initial values, we transform the (Li, Bi) coordinates of our reference stars into film coordinates ([^(x)]i, [^(y)]i). The distance si 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 si2. 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 aij and bij 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

The complete processing involves a number steps which are briefly explained here:
  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.

  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 9th or 10th magnitude. For Unix operating systems, XEphem is an excellent program which has support for the PPM catalog (limiting magnitude approx. 9m) and the Hubble Guide Star Catalog (limiting magnitude approx. 14m)
  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

The processing steps described in the previous section are shown in Fig. 5. The original image was a 28 mm wide-angle photograph of the summer Milky Way on Kodak Ektapress PJM film. Note how the rectangular image is transformed into a barrel shape.

example.jpg
Figure 5: Processing steps: (1) raw image, (2) background smoothed, (3) flat-fielded and contrast-stretched; the reference stars are marked as green dots, (4) transformed image.

5  The final image

Between 1997 and 2000, a total of 51 images were combined into a panorama covering the entire sky. The pictures were taken with a Minolta 28 mm lens riding piggyback on my Super Polaris DX mount. My observing sites were places as far apart as California's White Mountain Range and Cederberg Observatory in South Africa's Western Cape Province. Appendix A shows a complete list of the individual exposures.

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
In addition, articles on the image have been published in several magazines[3,4,5,6]

A  Exposure Table

No. Region f [mm] date UTC Film Site
1 Sco - Sgr - Sct 28 1997-07-06 07:31 - 08:00 Kodak PJM 2 1
2 CMa - Pup - Vel 28 1998-03-25 20:04 - 20:44 Kodak PJM 2 3
3 Vel - Cru - Car 28 1998-03-28 01:13 - 01:53 Kodak PJM 2 3
4 Sco - Cen 28 1998-03-28 02:25 - 03:00 Kodak PJM 2 3
5 CMa - Pup - Vel 28 1998-03-29 20:42 - 21:22 Kodak PJM 2 3
6 Sco - Sgr - Sct 28 1998-03-30 02:13 - 02:53 Kodak PJM 2 3
7 Sct - Lyr 28 1998-09-16 04:50 - 05:35 Kodak PJM 2 2
8 Cas - Aur 28 1998-09-16 08:27 - 09:12 Kodak PJM 2 2
9 Lyr - Cep - Cas 28 1998-09-17 04:59 - 05:44 Kodak PJM 2 2
10 Aql - Del - Aqr 28 1998-09-16 07:31 - 08:16 Kodak PJM 2 2
11 UMi - Aur 28 1998-09-16 09:19 - 11:00 Kodak PJM 2 2
12 Equ - Peg - And 28 1998-09-17 05:50 - 06:35 Kodak PJM 2 2
13 Aur - Ori - Mon 28 1998-09-17 09:34 - 10:19 Kodak PJM 2 2
14 Mon - CMa 28 1998-09-17 11:00 - 11:30 Kodak PJM 2 2
15 Aur - CMa 28 1998-09-19 11:42 - 12:22 Kodak PJM 2 2
16 Oph - Lyr 28 1998-09-16 05:44 - 06:29 Kodak PJM 2 2
17 Lyr - UMi 28 1998-09-16 06:37 - 07:22 Kodak PJM 2 2
18 And - Ari 28 1998-09-17 07:32 - 08:17 Kodak PJM 2 2
19 Aqr - Ari 28 1998-09-19 06:46 - 07:31 Kodak PJM 2 2
20 PsA - Cet 28 1998-09-19 07:38 - 08:23 Kodak PJM 2 2
21 Sgr - PsA 28 1998-09-19 03:37 - 04:22 Kodak PJM 2 2
22 Cet - Ori 28 1998-09-19 09:40 - 10:25 Kodak PJM 2 2
23 Cet - Lep 28 1998-09-19 10:31 - 11:16 Kodak PJM 2 2
24 Gem - Aur 28 1999-01-19 20:30 - 21:00 Kodak PJM 2 7
25 CMa - Pup - CMi 28 1999-01-19 21:30 - 21:56 Kodak PJM 2 7
26 Leo - Hya - Sex 28 1999-04-18 20:43 - 21:11 Kodak PJM 2 8
27 Leo - LMi 28 1999-04-18 21:17 - 21:45 Kodak PJM 2 8
28 Leo - LMi - UMa 28 1999-04-18 21:55 - 22:23 Kodak PJM 2 8
29 Leo - LMi - UMa 28 1999-04-19 21:31 - 22:04 Kodak PJM 2 8
30 Vir - Crv 28 1999-04-19 22:10 - 22:40 Kodak PJM 2 8
31 Leo - UMa 28 1999-04-19 22:59 - 23:27 Kodak PJM 2 8
32 Boo - Oph - Sco 28 1999-04-19 23:40 - 00:05 Kodak PJM 2 8
33 Sco - Sgr 28 1997-07-04 21:52 - 22:42 Kodak PJM 2 1
34 Sco - Sgr 28 1999-10-05 19:44 - 20:20 Kodak PJM 2 4
35 Sco - Cen 28 1999-10-06 18:06 - 18:47 Kodak PJM 2 4
36 Sco - Oph 28 1999-10-06 18:59 - 19:22 Kodak PJM 2 4
37 Pav - Tuc 28 1999-10-08 21:11 - 21:51 Kodak PJM 2 3
38 CMa - Car 28 1999-10-09 02:27 - 03:02 Kodak PJM 2 3
39 Pup - Vel - Cen 28 1999-10-10 02:30 - 02:57 Kodak PJM 2 3
40 PsA - Eri 28 1999-10-09 19:18 - 19:58 Kodak PJM 2 3
41 Cet - Eri 28 1999-10-09 20:07 - 20:47 Kodak PJM 2 3
42 Eri - Tuc - Car 28 1999-10-08 22:01 - 22:33 Kodak PJM 2 3
43 Mon - CMi - Cnc 28 2000-01-02 08:44 - 09:22 Kodak PJ 400 6
44 UMa 28 2000-01-02 10:10 - 10:50 Kodak PJ 400 6
45 Leo - Com - Boo 28 2000-04-02 21:14 - 21:59 Kodak PJ 400 3
46 Vir - Crv 28 2000-04-02 22:10 - 22:55 Kodak PJ 400 3
47 Boo - CrB 28 2000-04-02 23:04 - 23:49 Kodak PJ 400 3
48 Cen 28 2000-04-06 00:28 - 01:13 Kodak PJ 400 5
49 Sco - Boo 28 2000-04-06 01:30 - 02:15 Kodak PJ 400 5
50 UMa - Boo 28 2000-04-28 22:46 - 23:16 Kodak PJ 400 8
51 Oph 28 2000-04-28 23:20 - 23:45 Kodak PJ 400 8

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


File translated from TEX by TTH, version 2.72.
On 4 Feb 2001, 19:44.