Combining images

Sample mean as best estimator of parent population mean in least-squares, maximum likelihood sense. However, the mean is not especially robust in the case of outliers. Outliers occur in lots of astronomical examples, e.g., cosmic rays or filtering of stars. Combining images while doing the best job of rejecting outliers is a critical part of many data reduction/analysis tasks.

Some more robust estimators: median, but it produces larger error of the mean than does the mean. Often can make use of a priori knowledge about outliers: e.g., stars and CRs are always positive. This leads to routines like reject the maximum pixel. However, this leads to biases for all pixels without outliers. Hence, a better technique is min-max rejection; even this still leads to biases in pixels which have an outlier, and really throws away signal on others. Probably the best bet is to do n-sigma rejection, then recomputation of the mean. Problem here is that estimator of sigma can be very biased in the presence of outliers. So this really only is easily applied in the situation where error model is known. However, it can be applied more generally if you compute both mean and variance from the sample with the maximum value removed.

What if you don't have equal weight exposures, e.g. twilight flats? You need to normalize and potentially, consider the effect of different noise levels.

Note that the subtlety of combining frames comes into its full play when you are doing sky subtraction frames. Clearly these frames will have objects and these need to be filtered out. To do this, you must take many different sky frames with pointing shifts between them and combine the stack to make a sky frame. Be very wary of how you do this combination! It is very easy to do a straight median which to your eye will look superb; you will not see any residual sources. However, make sure to try the experiment of comparing several supersky frames by dividing them. Unless you are very careful in your combination, I'll place bets that when you divide your superskys that you will see remnants of wings of your objects. These will systematically corrupt your measurements of your program objects.

The method we've developed to optimize rejection of objects is the following. First, determine the shifts between all of the sky frames. Shift each of the frames to a common pointing and stack them up (we'll discuss how to shift images shortly). This provides a maximal S/N frame on which you can recognize faint sources. Identify all of the sources on this frame, probably automatically by some sort of positive peak detection. Make masks around each of the peaks, and make these masks larger than the peaks themselves so you also remove low level wings. Then shift each of these masks back to the pointings of each sky frame, and combine all of your sky frames at their original pointings, but ignore masked pixels in each frame. This method allows you to find and remove the very faintest objects you can see. If you were just to do peak rejection on each individual frame, you'll miss the faintest objects. If you just combine the frames without any masking at all, you'll get poor results. Note this technique is almost required if you find that the sky frames are changing with time. In this case, you'll just have to use the nearest neighbor sky frames (possibly just two frames!). However, you can find objects just fine using many more frames.

Measuring positions of objects and relative pointings between frames

It's clear that to do IR processing, you will be required to be shifting and coadding frames, because all of your object frames will not be at exactly the same pointing since you will be going off and observing sky and the telescope pointing won't repeat perfectly. In addition, you probably want to plan or making the pointings somewhat different for each observation so you can ignore bad pixels without losing pieces of your image. There are many other tasks for which you might need to shift or resample images, for example, if you want to combine images taken at two different telescopes, or two different wavelengths (maybe even with the same instrument if there are refractive elements!).

Before performing shifts, you need to determine how much to shift images. Generally, the quickest way to do this is to find some point sources on your image, measure their positions on each frame, and average the differences to get the frame to frame shifts. Beware if shifts are very large, because then you need to check to make sure distortion is not a problem; the presence of distortion leads to a different shift at different locations in the field. Aligning distorted images is much more complicated.

There are several ways to measure positions. A simple and accurate way to form the first moment of an image distribution:

$$<x>={\sum x I_{ij} \over \sum I_{ij}}$$

$$<y>={\sum y I_{ij} \over \sum I_{ij}}$$

Other possibilities are fits to the image shape. For example, fit a Gaussian to the image, solving for the brightness, width, and centroid. This requires a 2D non-linear least squares fit, which is somewhat more complicated, but still fairly simple. Other alternatives are Gaussian fits to the marginal distribution, a technique used the by the astrometry group at the US Naval Observatory, which has a large stake in determining accurate positions, since they are one of the main producers of quality astrometry in the world. Clearly, these methods only work if you have point sources with sufficient signal-to-noise.

An alternate way if you have extended objects or if you have lower S/N is to use cross correlation. Take your two images, multiply them together and sum the result. Do this multiple times with various shifts between the images. Fit a surface to the resulting sums as a function of row and column shifts, and find the peak. This is the desired shift. Of course, to do cross corrlations, you need to know how to shift images!