chi and h site logo

Home
About χ & h

Webcam amateur astronomy

Object
Detector
Optics
Camera
Mount
Computer and operating system
Driver software
Acquisition software
Data reduction
Presentation
Theory and tests

Data reduction

I use Starlink software for reduction as far as possible. Starlink is a project financed by the United Kingdom's Particle Physics and Astronomy Research Council (PPARC). The Starlink applications come in a number of packages. I use only a few applications from a few packages:

The software I add to this will often read and write Starlink NDF format, or it may use the libpgm or libppm libraries from Netpbm or the cfitsio library from NASA's High Energy Astrophysics Science Archive Research Center at the Goddard Space Flight Center. This may be an interesting point for Windows users, as it would be possible to port such programmes to Windows (by building them in a Cygwin environment).

I don't go into detail regarding the programmes I've written for these data reduction tasks. They are not in a releasable state, not named in a consistent manner, etc. Some day, may be.

During data reduction raw frames are grouped, aligned stacked and mosaicked. We should understand what I mean by these terms.

Frame
An individual image from the webcam.
Frame group
Some data aquisition applications like AstroVideo and Vega can add a number of frames - a frame group - immediately and without checking their relative alignment. The application will then write to disc not each frame but only the added frame group. Even if during aquisition individual frames are recorded to disc, during reduction one could blindly average frames into frame groups.
Stack
Frames or frame groups that supposedly show the same field of view, but actually do not due to tracking errors, lack of tracking or wind buffeting, are processed into a stack. This involves determining the shift - or shift and rotation - between frames or frame groups, and then to average them after applying the shift.
Mosaic
The field of view may be too small or it may be drifting too fast if there is no tracking. In that case several stacks need to be combined into a mosaic. The stacks may cover adjacent areas of sky with significant overlap, or they may try but fail to cover the same area of sky. Similar to stacking, the shift (and rotation) need to be determined and the stacks combined into a larger mosaic.
Capture
A capture is an entity of data recording, usually it comprises the frames that form a stack. But you could also imagine taking a long AVI recording during which you move the field of view to different areas of interest. Still, in the context of data reduction, a capture is the same as a stack.
Sub-capture
The CCDPACK software is limited to working on at most 100 frames or frame groups at a time. Hence it may be necessary to divide a capture or stack into several sub-captures of less than 100 frames. It may sound odd that professional astronomers put up with software that handles only 100 frames. But professional astronomers take very long guided exposures, they don't have to piece together their exposure from 100 or 1000 frames.
Target
A target capture is one taken on a sky object of interest. This is different from a bias capture, a dark capture or a flat field (cf. Theory / Detectors). I usually take a combined bias and dark capture with the same webcam settings as for the target.

Colour imaging with qastrocam

Renaming

The starting point of the reduction is a present working directory that contains one sub-directory for each capture. (The present working directory is the point in the directory tree that your command line terminal was located in with the "cd" command before you ran qastrocam.) qastrocam will have given it a name like brute-2003.02.01-10:14:04. Each capture contains one PPM file for each frame of the capture. The frame file names are 1.ppm, 2.ppm, ... 99.ppm, 100.ppm (here for 100 frames per capture).

The capture names are not bad, because their alphabetic order is a chronological order, and this can help you remember which capture was taken on which object or is a dark for which other capture. But I would say the first reduction step is to rename the captures to something more memorable, like "jup1" or "dark5". Of course, now you have to remember for the new names which dark belongs to which targets.

mv brute-2003.02.01-10:14:04 moon01
mv brute-2003.02.01-10:14:34 moon02
mv brute-2003.02.01-10:17:10 dark1
mv brute-2003.02.01-10:14:04 jup1
mv brute-2003.02.01-10:17:10 dark2
mv brute-2003.02.01-10:14:04 sat1
mv brute-2003.02.01-10:17:10 dark3

My second reduction step is to rename the frame names from plain numbers to zero-padded numbers. Assuming that no capture has more than 9999 frames, this should change the names from 1.ppm, 2.ppm, ..., 99.ppm, 100.ppm to 0001.ppm, 0002.ppm, ..., 0099.ppm, 0100.ppm:

for dir in moon01 moon02 dark1 jup1 dark2 sat1 dark3 ; do
  cd $dir
  for file in ?.p[gp]m ; do mv $file 000$file 2>/dev/null ; done
  for file in ??.p[gp]m ; do mv $file 00$file 2>/dev/null ; done
  for file in ???.p[gp]m ; do mv $file 0$file 2>/dev/null ; done
  cd ..
done

This is for a Bourne-type shell like sh or bash. For a C-type shell like csh or tcsh the for ... ; do ... ; done construct has to be written differently and the output re-direction may also have to be different.

Weeding

Now comes the really tedious bit, frame weeding. We want to remove frames that are not good. They might suffer from cloud passing through or the wind may have been shaking the telescope to blur the image. But most of all, the image quality at long focal lengths comes and goes: There are some wonderful frames, a number of good ones and a majority of mediocre or bad ones, where the amount of detail is just unsatisfactory. AstroVideo can find the best frames in an AVI automatically, but usually we have to use our own eye and judgment. A reasonably convenient way of doing this for qastrocam frames is to start a background process for displaying all frames with ImageMagick (comes with Red Hat Linux) and then to run a query-and-remove command to delete and keep the frames as we want:

/usr/bin/X11/display *.p[pg]m &
rm -i *.p[pg]m

If you forget the -i all your frames will be gone without question!

Master darks

The first data to process are the dark captures, because they do not depend on any other captures. Second would be the flat captures, if I took any. Last are the target captures, which need dark subtraction (and division by flat, if available).

I've written a programme makedark to process the frames of a dark capture into a master dark cube. It uses the libppm graphics library from Netpbm to read the frames direct from the PPM files and it uses the cfistio library to write out the master dark in FITS format. In this process the data are re-sorted so that the output cube contains three images for RGB. The output file is called dark.fit.

makedark looks for frames 0001 to 1000 and uses those it finds. They are simply averaged pixel by pixel. While the input numbers are integers from the set {0, 1, ..., 255}, the output numbers are floating point numbers in the interval [0,255]. Consequently the output FITS is BITPIX=-32, meaning IEEE 4-byte floating point numbers.

for dir in dark1 dark2 dark3 ; do
  cd $dir
  makedark
  cd ..
done

This, in my view, is an important point. Astro webcam'ers who work under Windows tend to have a fixation on integers and on 8-bit integers for that matter. True, the input is 8-bit digitised and the final image you show to your friends is also a graphics with 8-bit integers. But it is important to do all the reduction in floating point numbers and with 16 or 32 bits so that you don't loose in the arithmetic what you were trying to gain by taking many frames. (Cf. Theory / Calibration)

Stacking

Similar to makedark, I've written the programme stackdisc. This looks for frames 0001 to 1000 of a target capture, aligns and averages those it finds and writes out the result as floating point FITS. It also requires a FITS cube called dark.fit. Symbolically link the appropriate master dark before running stackdisc:

cd moon01
  ln -s ../dark1/dark.fit dark.fit
  stackdisc
cd ../moon02
  ln -s ../dark1/dark.fit dark.fit
  stackdisc
cd ../jup1
  ln -s ../dark2/dark.fit dark.fit
  stackdisc
cd ../sat1
  ln -s ../dark3/dark.fit dark.fit
  stackdisc
cd ..

The output file is called stack.fit and covers the same field of view as the first frame found.

Why do I call it stackdisc and not just stack? The automatic alignment of frames is not so easy to do right. stackdisc does not work well with star fields. It is written for the case where we have a bright disc entirely in the field of view. The disc does not have to be round, so Jupiter or Saturn will be fine. stackdisc also happens to work well where the disc is larger than the field of view, such as a detail on the Moon or sun spots.

But stackdisc can fail if the limb of the Moon or Sun run through the field of view. When this becomes a problem, I will write a similar stacklimb programme. In fact I already have a stackstar (cf. Reduction / Grey / Stacking)

Sun and Moon: monochrome mosaic

For the Sun and the Moon colour is not important. The simplest way to proceed is to pick one of the three colours and keep that as a monochrome image. To prepare for the further reduction this is converted from FITS to the Starlink format NDF:

for dir in moon01 moon02 ; do
  cd $dir
  fits2ndf stack.fit stack
  subset "stack(,,2)" grey
  rm stack.sdf
  cd ..
done

One might be tempted to average RGB into grey, just to use more of the raw data. But we then suffer from differential refraction. At long focal lengths the B image is visibly higher in the sky than the R image: the Earth's atmosphere acts like a prism and splits the light into a spectrum of its component colours.

For the Moon it is likely that I have a number of stacks that show different parts, but pairs of stacks will overlap with each other. Say, I might have a sequence of stacks to cover much of the terminator. The intention is then to make a mosaic of these stacks.

I need to form groups of stacks such that within the group there is one stack that has sufficient overlap with all other stacks in the group. I can then use my regimag programme through its cross2.sh wrapper script. This will display the reference stack and one of the other stacks and I indicate with the mouse the approximate alignment. cross2.sh will refine the alignment and store it with the stacks. cross2.sh goes through all NDF images in the present working directory. After cross2.sh I run the makemos application from CCDPACK, which will construct a larger image and insert and average the shifted constituent stacks.

mkdir moon101
cp moon01/grey.sdf moon101/moon01.sdf
cp moon02/grey.sdf moon101/moon02.sdf
cd moon101
  cross2.sh moon01
  makemos in="'moon??'" out=grey scale=y zero=y method=mean
cd ..

Inspect the result. It is possible that the scale/zero fit is not a good idea. In that case re-run makemos with scale=n zero=n.

If there are many stacks to combine, you have to first build up local mosaics from stacks that all overlap with their reference stack, then combine neighbouring local mosaics into regional mosaics, and so on until you arrive at the global mosaic.

Planets: colour images

For planets it is likely we want to keep the colour. In that case I extract all three colours as NDF images and align them with each other. This takes care of differential refraction. For the alignment I have a programme that works on NDF images. Before running its align.sh wrapper we have to remove the file stack.sdf, because align.sh will align all .sdf files it finds, and it must find only the three RGB images.

for dir in jup1 sat1 ; do
  cd $dir
  fits2ndf stack.fit stack
  subset "stack(,,1)" stackR
  subset "stack(,,2)" stackG
  subset "stack(,,3)" stackB
  rm stack.sdf
  align.sh stackG
  cd ..
done

align and align.sh do not stack, they only determine the shifts required to align the images and store this with the images. Before we do combine RGB into colour or grey, we should visually inspect one of the images and decide on a suitable smaller part to keep. We don't need the full 640x480 pixels, only enough to see the planet and a bit of black sky. Having decided on a pixel range, we can then extract that for each of the colours. A grey image can then be constructed as an average "mosaic", making use of CCDPACK's makemos.

To construct a 24-bit colour image I have written the truec programme. It takes three congruent NDF images and applies the same linear stretch to each. Don't stretch it too much: If you choose as upper limit the actual maximum of the input images the brightest part of the colour image will be white, i.e. show no colour. Better leave 10 per cent of the dynamic range unused. truec puts the image as ASCII PPM to standard output. Redirect it into a format converter to make better use of disc space.

cd sat1
  subset "stackR(101:340,201:450)" subR
  subset "stackG(101:340,201:450)" subG
  subset "stackB(101:340,201:450)" subB
  makemos in="'sub[RGB]'" \
    out=grey scale=n zero=n method=mean
  truec black=0 white=150 red=subR green=subG blue=subB | \
    pnmtopng > colour.png
cd ..

This fuller script loops through captures and colours to

  1. symbolically link the dark,
  2. subtract dark, align and stack,
  3. extract RGB and align them with each other,
  4. cut off the edge,
  5. unsharp mask (cf. Presentation / Unsharp masking),
  6. combine RGB with linear stretch.
for dir in jup1 sat1 ; do
  cd $dir
  ln -s ../dark1/dark.fit dark.fit
  stackdisc
  fits2ndf stack.fit stack
  subset "stack(,,1)" stackR
  subset "stack(,,2)" stackG
  subset "stack(,,3)" stackB
  rm stack.sdf
  align.sh stackG
  for col in R G B ; do
    subset "stack$col(6:635,6:475)" sub$col
    unsmsk2.sh sub$col mask$col 10 0.7
  done
  truec black=0 white=60 red=maskR green=maskG blue=maskB | \
    pnmtopng > mask.png
  cd ..
done

This does not cut down the image to the region of interest. This is done with the intent of combining a number of consecutive captures into a movie. One would copy all captures' sub[RGB].sdf into one directory, align the captures by determining the shifts between their G channels and then apply the same subset to all. This is what is done below. Notice how the first loop can work on directory name expansion from wildcards, while the second loop has to use actual file name fragments.

for dir in jup?? ; do
  for col in R G B ; do
    cp ${dir}/stack${col}.sdf jupall/${dir}${col}.sdf
  done
done

cd jupall
aligncol jup01G

for dir in jup01 jup02 jup03 jup04 jup05 jup06 jup07 \
    jup10 jup11 jup12 jup13 jup14 jup15 jup16 jup17 jup18 jup19 ; do
  for col in R G B ; do
    subset "${dir}${col}(200~250,232~250)" ${dir}sub${col}
    unsmsk2.sh ${dir}sub${col} ${dir}mask${col} 10 0.7
  done
  truec black=0 white=60 red=${dir}maskR green=${dir}maskG blue=${dir}maskB | \
    pnmtopng > ${dir}.png
done

cd ..

Grey imaging with qastrocam

Renaming and frame uniqueness

See also the corresponding section for RGB frames. The starting point is similar to the colour case. The frames are PGM instead of PPM files. Another difference stems from the different webcam. The QuickCam VC will have been set to maximum exposure time of 1.83 s, but qastrocam will have recorded a PGM frame every 0.6 s. Hence each actual frame has been recorded on average 3.05 times. We should check the recorded fames and discard such duplicates. If we do this before changing to zero-padded frame numbers this can be done with

for dir in dark1 jup1 jup2 jup3 sirius ; do
  cd $dir
  n=1
  while [ -r ${n}.pgm ]; do
    m=$(($n+1))
    if [ -r ${m}.pgm ]; then
      if diff --brief ${n}.pgm ${m}.pgm ; then
        rm -f ${n}.pgm
      fi
    fi
    n=$m
  done
  cd ..
done

The diff command returns an error code that tells whether the current frame and the next frame are identical or not. If they are, the current frame is removed. So only the last recording of each genuine frame remains after this.

Weeding

See also the corresponding section for RGB frames. For star fields at short focal length weeding is usually not necessary. If part of the capture is affected by passing cloud or someone shining a light into your lens then those frames should be weeded out, of course.

Master darks

See also the corresponding section for RGB frames. I have a second programme, makedarg to work from PGM frames and generate a dark image from them. Note that this one ends with "g" for grey (instead of "k" for kolour; no: kolour is not a proper German word).

Stacking

See also the corresponding section for RGB frames. I have a second programme, stackstar to align and stack grey images. The alignment algorithm is different, more suited to a star field, including one that drifts across the detector.

Mosaic

See also the corresponding section for Sun and Moon. If we have a number of stacks that overlap and jointly cover a larger area, we can form a mosaic. One way is to use the cross2.sh wrapper around regimag as we do above for Sun or Moon.

For star fields there is an alternative, as the CCDPACK package is written for this sort of thing. Unlike cross2.sh, this method does not require each stack to overlap with a single reference stack. So we can mosaic a larger area in one go. There is, however, a limit of 100 stacks that can be processed at once in this way.

The first task is to find objects in each stack. Have a glance over the output to make sure at least 3 objects are found in each stack, and also that not too many objects are taken, which may be insignificant or even false (noise peaks). Reduce the percentile to find more objects and vice versa.

findobj in='*' outlist='*.find' \
  autothresh=y override=y useper=y percentile=99.9 minpix=3

The second task is to have a first stab at aligning the stacks with each other. This is done by comparing the objects found.

findoff inlist='*' outlist='*.off' \
  error=3 minmatch=1 override=y accept

The third task is to determine the final pixel transform between the stacks. This is more than a simple shift by so many pixels in x and y, the shift includes fractions of a pixel and also a rotation. There is usually a slow rotation of the field when the stacks have been taken over half an hour or more.

register inlist='*' fittype=2 refpos=1 accept

Now each stack is transformed to a common grid of new pixels and then these resampled stacks are then averaged into the mosaic. The scale/zero fit in makemos takes time, and may also deteriorate rather than improve the mosaic. On the other hand, fitted scales far from 1.0 and large fitted zero level differences may indicate that something went wrong in the alignment.

tranndf in='*' out='*_t' method=linint conserve=y
makemos in='*_t' out=mosaic scale=y zero=y method=mean

Copyright © 2003 Horst Meyerdierks
$Id: reduce.shtml,v 3.3 2004/02/21 18:13:39 hme Exp $