Data Reduction: details and subtleties

Data reduction steps

Generally, remove instrument signatures in reverse order from which they were introduced. First consider basic processing steps applicable to all array detectors, i.e. regardless of whether instrument is imaging or spectroscopy. Possible correction steps include:

1. A/D correction: rarely done.

2. Bias level subtraction: always done for CCDs, usually not for IR detectors

3. Bias structure subtraction (superbias): often done, but not always warranted.

4. Linearity correction: rarely done with CCDs, typically done with IR arrays.

5. Dark current subtraction (superdark): occasionally done with CCDs, depending on detector, exposure times; usually done for IR arrays.

6. Preflash subtraction, or deferred charge correction: rarely done with modern CCDs.

7. Shutter shading correction: required if effect is significant.

8. Flat field correction: always done.

9. additonal complicated things: scattered light (has different flat field characteristics). Amplifier hysteresis correction (probably only good qualitatively!). Similarly, CTE correction. Hot pixel fixup. Residual image subtraction. Pixel area correction.

10. Fringing subtraction/sky subtraction: occaisionally done where warranted.

For each of these steps, it is important to understand whether the correction is additive or multiplicative.

Bias level subtraction

Use overscan/underscan if available. If not, you must use bias frames or dark frames of the same exposure time. Note IR detectors should be already bias-subtracted if they are using doubly-correlated reads, but often there is residual bias which is often exposure-time dependent, and if so, can be subtracted out with darks of the same exposure time (at the same time doing dark subtraction!). Warnings about overscan regions: amplifier hysteresis. Vertical and horizontal overscan. (vertical often not taken, rarely applied). Note possible drift of bias level with time (apparent in vertical overscan), so possibly fit. If fitting for varying bias, beware of bias jumps.

Bias structure subtraction

Readout electronics sometimes introduce fixed pattern noise, and in some cases, it is possible to subtract this from your data. However, you can only do this if the pattern is repeatable. You must check; beware of a ``fixed'' pattern which changes with, e.g., telescope pointing. Often, this noise is variable in time or location, and consequently you only make things worse by attempting to correct for it.

Check bias frame histograms. Should be zero mean and Gaussian with width given by readout noise. This gives you first observed estimate of readout noise. If it's not zero-mean or has significantly larger width, look for fixed pattern. Construct superbias frame from lots of bias frames. Repeat several times and see if histogram is repeatable.

Required number of bias frames depends to some level on how bad the bias structure is. If it's really bad, i.e., much worse than the readout noise, you'll improve things with only a few bias frames. However, usually, the bias structure is subtle and less than the readout noise. In this case, you can often add noise by doing a superbias subtraction. In the absence of any bias structure, the superbias will have noise of rn/$\sqrt{{N}}$ , and consequently, your superbias subtracted frame will have effective readout noise of $\sqrt{{rn^2 (1 + 1/N)}}$ . So probably want at least 10 frames, more if you really want to go faint and will be readout noise limited.

IR data should not required bias structure subtraction with doubly-correlated reads, but note that some instruments do!

Note that it is possible to do bias subtraction simultaneous with dark subtraction if you take calibration of the same length. This is typcially what is done with IR detectors.

Linearity correction

For some detectors, especially in the infrared, you must make a correction for the nonlinearity of the devices. To do this, you must measure the output intensity at a range of known input intensities. To do this, you must be able to accurately control the input intensities. This is probably best done by taking data at different exposure times; since shutter timing is usually quite precise, the ratio of the exposure times give the ratio of the input intensities. However, for this to be true, you must use a light source which is stable, i.e., does not vary with time (at the one percent level); this is often hard to come by.

To make a linearity curve, you take a series of frames of increasing exposure level, usually by increasing the exposure time. To do linearity, you must know the relative levels of your different frames, which is the ratio of exposure times if the light source is stable. However, most light sources are not stable. To do a linearity correction, you need at least a light source which does not vary on short time scales. Then you can intersperse lamp monitoring exposures to try to calibrate out the light variations. Generally you can do this by interspersing some fiducial exposure between each of your successive light levels. For example, a linearity sequence might consist of an exposure sequence like: 1s, 2s, 1s, 4s, 1s, 8s, 1s, 16s, etc. The longer exposures are used to make the linearity curve, and the 1s exposures are there to test for source variability; if the latter occurs, the short exposures can be used to correct the other exposures for changing source intensity (as long as the changes are of sufficiently long timescale that they are well sampled by the short exposures). When you are done, make a plot of light level in each of the fiducial exposures. Look for any drifts in the lamp brightness, and if they exist, you will need to calibrate them out of your linearity exposures. Plot the linearity exposures as DN/s vs mean count level. You need to do this in regions which has a roughly constant sensitivity; don't average over the whole array; it is even possible to do it on a per pixel basis. Fit some sort of curve (low order polynomial) through the data. Then for every observation, use this curve to correct the observed count rate to the count rate you would have received at some fiducial level.

At many observatories, the linearity correction may have been previously measured and will be available for you to use without needing to measure it yourself.

Note that the correction will be larger the further your program observation is from the fiducial light level. To minimize errors from errors in your linearity correction, take your program exposures as close as possible to a single light level; when exposure levels differ a lot, you will need to be trusting your linearity correction more.

Dark current subtraction

Detectors generate dark current arising from thermal excitation. Generally, the dark current in a device is different for each different pixel. Most pixels have very low dark current (because an operating temperature is chosen to bring dark current to an acceptably low level.) However, some pixels will have larger dark current (warm pixels), and some pixels may have very high dark current (hot pixels). The amount of dark current on a frame will obviously be larger for longer exposure times.

If the dark rate in a given pixel is constant in time, one can subtract the dark current by using dark frame, i.e. frames with the shutter closed but of finite exposure time (not bias frames), and making ``superdark'' frames by combining multiple dark frames (to beat down photon and readout noise). Sometimes the dark rate in a pixel does not appear to scale linearly with time, so sometimes you have to beware of using dark frames at some exposure time to correct object frames of a different exposure time; in the ideal world, you would take dark frames at each different exposure time you will be using for your object frames. However, dark frames are time-consuming to take. Also, you need to be careful about taking dark frames during the daytime, because sometimes instruments have small light leaks which make daytime darks useless. Cloudy nights are an ideal time to take dark frames!

For pixels with very low dark rates, which is often the majority of the pixels, you need to take very large numbers of dark frames to avoid being dominated by readout noise. Often, however, you only have time to take a relatively small number of darks, since they must be of long exposure time. Conseqeuntly, you might wish to consider a strategy where you make a superdark frame, but set all pixels for which dark rate is too small to be accurately measured to be identically zero. This way you avoid just adding readout noise to these pixels, but you still use the measured dark rate to correct for warmer pixels.

Preflash subtraction, or deferred charge correction

For CCDs with charge transfer problems at low light levels (mostly older devices), nonlinearities can be avoided by preflashing each exposure with a small pedestal of light, although this is not ideal as it adds background noise to your object frames. If you add a preflash, it must be subtracted during data reduction; to do so, one would take a series of preflash-only frames, and combine them to make a preflash subtraction frame.

If your object frames have sufficient background by themselves so that all pixels are above the nonlinear low levels, it is possible to correct for deferred charge effects without the use of a preflash frame. Details of how to obtain the required calibration data and construct the relevant calibration frames are given in Gilliland in ``Astronomical CCD Observing and Reduction Techniques'' (PASP Conf. Series 23).

Shutter shading correction

Since shutters take a finite amount of time to open, the true exposure time will be a function of location on the frame. Generally, shutters move quite fast, so this is only a measureable effect for short exposure times, and even then, only for detectors which are physically large (i.e. have large shutters). However, this is becoming more common with very large format CCDs. If you care about photometry at the one percent level on short exposures with such a device (and you often do, because standard star exposures are often relatively short), you may need to consider a shutter shading correction. This can be derived by comparing N short exposures with exposure time t with a single long exposures of exposure time Nt ; for example 20 1s exposures with 1 20s exposure.

For any frame, the observed number of counts, O , will be given by:

O = S(t + $$\Delta$$t(x, y))

where S is the true count rate, t is the commanded exposure time, and $\Delta$t gives the difference in exposure time at different positions on the detector. From this you would get:

NO_t = NS(t + $$\Delta$$t)

O_Nt = S(Nt + $$\Delta$$t)

If a sufficiently long exposure is taken in the latter case, then, to first order,

O_Nt $$\approx$$ SNt

so

$${N O_t \over O_{Nt}}$$ $$\approx$$ $${(t + \Delta t)\over t}$$

so from a ratio of the two frames, one could compute $\Delta$t , and use the first equation above to correct every observed frame to the true count rate for the commanded exposure time:

S = $${O t\over t+\Delta t}$$

Flat fielding

Flat fielding is the calibration step which corrects for varying throughput as a function of position in the field. This arises from variations in pixel-to-pixel sensitivity, illumination patterns, and filter throughput variations, i.e. it is not necessarily just a detector effect. All of these source cause multiplicative field-dependent errors.

In principle, these are easy to correct for; all one needs to do is to observe a source which is spatially flat, and any deviations from a flat response give the field-dependent response variations. One then corrects all program observations by dividing by this flat-field frame. However, for accurate correction, one needs to observe a perfectly flat source, which can be difficult. Often, a white spot on the inside of the dome is used which is way out of focus, but, at some level this provides a different illumination than sky because it is a near-field target. Also,one needs to consider the color dependence of the flat field response; you need separate flats for every filter. For broad filters, you also need to consider the color of the light source and how it compares with the color of your objects.

Individual pixel response variations require high S/N to measure; generally, you should aim for accuracy of a small fraction of a percent. However, be aware that computed S/N doesn't give true accuracy for large spatial scales, because the intrinsic limitation is the uniformity and color of light source. The main problem in flat fielding is systematic errors, not random errors.

Generally, flats are constructed using one of three techniques, or some combination of these three: dome flats, twilight flats, and dark sky flats.

Dome flats are good for measuring the pixel-to-pixel variation if color match is good or if your array doesn't have strong color dependence; often, however, color match is poor. Because of illumination problems, however, dome flats are usually limited to a few percent accuracy for large scales.

Twilight flats can be good for pixel-to-pixel if you don't have too many filters and if they're not too narrow and if it doesn't take too long to read array. They are generally much better than dome flats for large scale structure. They can have some color mismatch problems. The biggest problem with twilight flats is that twilight is short, so it can be difficult to get data in all filters; certainly, one usually needs to use both evening and morning twilight. A lesser problem is that one will occasionally see stars in the twilight flats; this requires that you take several frames in each filter, moving the telescope between frames, and then reject stars in the process of combining the multiple frames. Additional problem is if you miss twilight, due to weather or scheduling.

Alternatively you can use your dark sky exposures to make the flats. However, generally, the level level is low (in the optical) so you can't build up high S/N. Consquently, when using these ``superflats'', one determines the pixel-to-pixel response with dome flats, and then you get the large scale variations with "superflats". Combine many frames from the night rejecting stars and fit surface to get large scale variation. A great advantage of this technique is that you match the color of the night sky perfectly.

What are the effects of systematic flat fielding errors? For no backgrounds, the fractional error in your measurement is the fractional error in your flat. However, in the real world, we need to subtract off background, and the presence of flat-fielding errors make you estimate the sky wrong, if you measure the sky at a different location from your object. You observe S + B. If you subtract xB, this gives S + (1-x)B, which differs from S by:

error = $${S + (1-x)B\over S}$$ = 1 + (1 - x)*$${B\over S}$$

so, for example, if S = B and x = 0.99 (one percent error in sky), then you make a one percent error. However, if B = 100S , then you make a factor of two error!

Thus flat-fielding becomes critically important if you are measuring sources which are fainter than the sky level. This occurs in the optical as you go very faint, and almost always occurs in the IR, because the background is so great.

Note that the same problem occurs if the background is intrinsically spatially non-uniform. This can be the case in the IR because of thermal contributions from the telescope/dome itself. This is especially troublesome because the background non-uniformity may very well change with time and telescope pointing.

What can you do about this? If all your frames are subject to the same flat-fielding or additive background errors, then you can beat the problem by observing the sky on the same pixels as your object. If you then subtract the sky frame from the object frame, you end up with a difference frame which has a much lower sky level as compared with your object (but still possibly non-zero if the background level varies in time). Consequently, the errors introduced by the flat-fielding errors go down. In fact, if the sky subtraction is perfect, you can reduce your object errors down to the percentage error in the flat field. However, remember that in doing this you introduce noise from photon statistics, since the subtracted frame has $\sqrt{{2}}$ the error of the original frame (if the sky dominates). So you don't necessarily win by doing this every time: only if the sky is brighter than the object. Depending on your requirements, you might be able to get away with observing the sky for somewhat less time than the object, if you want to keep the big systematic errors small and don't care about losing S/N at the faintest levels; but beware, you may pay something for your time saved. You can win this back if you average several sky frames for each object frame, but you have to beware for the sky frames changing with time.

Note if your object is small, you can do this without losing any observing time by dithering the object around the field. If your object is large, however, you increase your overhead substantially. But you don't have any choice.

Flats in the IR

Flats are additionally complicated in the infrared, because every frame you take, including flat field frames, have an additive extra component arising from the radiation coming from the telescope and dome structures. This is more of a problem at longer IR wavelengths, because at shorter ones, the sky dominates. Given this, how do you make a flat at all? One possibility is dome flats with and without lights: subtract the two for your flat. Of course, this has the potential problems associated with dome flats: incorrect illumination and color mismatches. But what can you do? Take both this kind and sky flats. If the variable spatial additive component from your telescope is small in the sky flats, you'll probably do better with these.

Since we are talking systematic errors in sky determination causing the biggest problems in your photometric accuracy, it is very difficult to estimate the errors analytically. Remeber, systematic errors are the most troublesome! Pretty much the only way you can check the accuracy of your flats is to take images of the same object at several different field locations and empirically determine how well the measurements repeat. This is highly advised in the IR, at least once with every instrument.

Scattered light and flat fielding

Another problem sometimes encountered with flat fields occurs if there is light that makes it to the focal plane without having come through the normal optical path; this is generically called scattered light.

If there is scattered light in object frames, it may be possible to remove it as an extra background contribution (that is not uniform across the frame).

But if there is scattered light in your flat fields, it can cause significant problems for photometry in all of your frames since you'll divide your frames by an incorrect field-dependent response function. Scattered light in flats can be espeically problematic since flats are often taking when there is significant light in the dome.

You may be able to check for the presence of scattered light by looking at flats taken with the telescope in different positions: in particular, if you can take flats at several different position angles, comparing these often can indicate whether there is a problem.

Flat fielding and distortion: pixel area effects

In many modern instruments, the field of view is sufficiently large that optical distortion can be significant, leading to a different pixel scale at different locations in the focal plane. This can have a significant effect on flat fields and requires you to carefully consider your flats.

If different pixels subtend different areas on the sky, the pixel area changes will be reflected in your flat field, so if you divide by a raw flat field, you will also be normalizing out different pixel area effects. But if you are observing unresolved sources, the light will just be spread out over a different number of pixels (i.e. a field dependent point-spread function in pixel coordinates), so if you divide by a flat field with the pixel area effects in it, you'll erroneously change the relative brightnesses of stars across the field-of-view!

If this is a significant effect, you will want to either remove the pixel area contribution from the flat-field (which you can do if a distortion solution has been obtained for the instrument), or you can resample the data onto a grid ofequal pixel areas (but this has the disadvantage of introducing correlated noise).

Drift-scanning and flat fielding

One technique for reducing errors from flat fielding: drift scanning (also called TDI, time-delay integration mod). Mechanics of drift scanning. Implications for improving flat fielding. SDSS as an example. Time variation along the strip: issues/opportunities for photometry and moving objects.

Constructing calibration frames : MOVED TO LAB

Combining images

Many astronomical applications require one to combine a stack of images. Several examples of this are for data reduction, e.g., the construction of flat fields, superdarks, superbiases, and sky frames for IR observing. Additionally, one may want to combine object frames, e.g. in IR observing when many frames are taken at different dither positions. For the data reduction products, multiple frames are taken to improve S/N and also to reject outliers in a single observation, e.g. cosmic rays and/or stars.

The best estimator of parent population mean in a least-squares, maximum likelihood sense, is the sample mean. However, the sample mean is not especially robust in the case of outliers. Outliers occur in lots of astronomical contexts, e.g., cosmic rays, filtering of stars, or just bad data. 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 (about 25% larger for normal distribution) 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 maximum-pixel rejection. 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; may work better if you compute both mean and variance from the sample with the maximum value removed). Alternatively, apply using error model; for example, compute sigma from the median value and a noise model, use this to reject outliers, then average the remaining data points .

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

Biases

Take multiple bias frames: Number depends on nature/amplitude of fixed pattern noise, readout noise. Can be done during daylight.

Combine frames to make ``superbias''. Each frame should be bias (overscan) subtracted first. Combine frames using technique to reject outlying pixels, e.g. median. (IRAF: imcombine, zerocombine)

Inspect superbias for fixed pattern noise, any net signal after overscan subtraction. Compare suberbiases taken at different times/positions.

Darks

Number depends on nature/amplitude of dark current. Need to insure dark conditions to get necessary data. Use cloudy nights if available. Long exposures necessary, especially if dark current is low.

Combine multiple frames to make ``superdark''. If dark current is low, consider clipping low values (lower than expected readout noise) to zero to avoid adding noise. (IRAF: imcombine, darkcombine).

Inspect superdark to determine whether dark subtraction is warranted. Compare individual darks with superdark to see how repeatable dark current is. If possible, compare dark frames taken at different exposure times to see how well linear scaling with time works.

Flats

Take multiple methods (dome + sky) if at all possible. Remember to dither telescope for on-sky flats. Take at different positions, position angles if possible. If taking on-sky, dither between exposures.

Combine multiple frames in individual filters. Normalize frames to common value first (esp. important for twilight, night sky frames!). Normalize final flat to unity to preserve noise characteristics after flat fielding. (IRAF: imcombine, flatcombine).

Compare individual flats with superflat to check for flat stability. Compare flats taken at different times/positions.

If significant distortion is present, consider removing pixel area contribution from flat field.

Other items

Checking detector behavior. Image histograms

Quick check of data quality from image histograms.

Linearity

Usually IR detectors only. Take linearity sequence, as previously discussed. Plot mean count rate as a function of number of counts. Fit some function to this, to be used to correct frames. Note that if illumination varies significantly across the field, may wish to split field into multiple sections.

Determining readout noise and gain, linearity check

Consider an ensemble of measurements taken at a light level L. The noise in this ensemble should be $\sigma^{2}_{}$ = LG + rn² , where G is the gain and rn is the readout noise in electrons, and $\sigma$ is the noise in electrons.

If you do this at a lot of different light levels, then you can plot $\sigma^{2}_{}$ vs L , and the slope should give you G and the intercept rn² . However, remember that if you compute $\sigma$ from the images, this gives $\sigma$ in DN, so the slope will give you 1/G . This test is also excellent for checking the basic performancer of a detector. Deviations from nonlinearity can also usually be seen on such a plot.

However, in its most straightforward application the test is very time consuming and hard to analyze: you have to take many exposures at each different light level, and then determine a gain and readout noise for each pixel and look at them all. It is much easier just to use the set of all pixels as your ensemble at each light level. However, you can't do this directly, because each pixel may have a different sensitivity and different fixed pattern noise, so you're not measuring a true ensemble. If there is significant variation of sensitivity than you can't use a whole area at all, because the noise properties will vary across the area. You can avoid these problems by working with the difference between pairs of observations: if the light level is the same in the two images, then you'll be left with an image that only has noise.

Specifically, take a pair of images and form the difference. The expected noise is

$$\sigma^{2}_{}$$ = 2(LG + rn²)

where L is the light level, G is the gain, and $\sigma$ is the noise in electrons. Since $\sigma$(electrons) = G$\sigma_{{DN}}^{}$ , we have

$$\sigma_{{DN}}^{2}$$ = 2($${L\over G}$$ + $${rn^2\over G^2}$$)

. You can directly measure $\sigma$(DN) from your difference image. Make sure to do it over a region which doesn't vary significantly in light level. Now repeat the measurement at a variety of light levels and plot $\sigma^{2}_{}$(DN) vs L . If you fit a line through this, the slope is 2/G and the intercept is 2rn²/G² . If a straight line doesn't fit the points, then there is some sort of problem, which you should probably track down.

You can abbreviate this test if you just want to get an estimate of the gain and readout noise. First, take a pair of bias frames. These have light level of zero, so the noise from the difference just gives you $\sqrt{{2}}$rn . (Note that you still need to take a pair in case there is superbias structure). Then take a pair at a high light level; at this level the readout noise is probably negligible, and you can determine the gain from

G = 2$${L\over \sigma_{DN}^2}$$

Basic data reduction techniques: MOVED TO LAB

Optical observing and basic reduction.

• subtract bias (overscan)
• subtract superbias
• subtract dark current (scaled to appropriate exposure time)
• shutter shading correction (scaled to appropriate exposure time)
• flat field

IR observing and basic reduction.

• bias/dark subtract: darks should be taken at each different exposure time that you use for your program objects
• sky subtract
• flat field

Making IR sky frames, where one wants to combine frames where you need to filter out astronomical objects (e.g. stars) can actually be rather tricky, as the profile of objects include pixels of a wide range of intensities, going all the way down to zero intensity; the fainter pixels in the wings of profiles are not dramatic outliers when one considers a single pixel at a time! All of the combining techniques discussed above consider a single pixel at a time, and don't take advantage of knowledge that outliers in some pixels will be correlated with outliers in adjacent pixels. 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, if you try the experiment of comparing several supersky frames by dividing them, you will likely see remnants of wings of your objects when you divide your superskys. These will systematically corrupt your measurements of your program objects.

To filter out stars, etc., the best bet is probably to identify the object, then mask out large regions around each object before combining; you then use your knowledge of the location of objects and the profile of each object rather than making a decision based on individual pixel values.

The same issue of not being able to identify marginal outliers in the wings of stars applies to the ability of being able to identify the faintest stars. For optimal results, consider combining all of the frames to detect the faintest objects, and use the combined frames to construct your masks:

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.

Coadding images

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 to make the pointings somewhat different for each observation so you can reject bad pixels without losing pieces of your image. There are many other tasks for which you might want to shift or resample images, e.g. combining multiple observations of the same object.

Coadding images once they are aligned is just simple image arithmetic, so the new image processing task requires is the ability to register images, which has two parts: measuring the shifts between images, and then shifting them to align them.

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 more complicated; the same goes for images that are rotated relative to each other.

Measuring positions of objects and relative pointings between frames

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.

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. Often, the shifts for your cross correlation can just be simple integral pixel shifts; you can derive fractional shifts by fitting for the peak of the cross-corrleation function.

Cross correlation is also widely used in spectroscopy to measure velocities. One needs to remember that a shift from velocity is not a constant shift in wavelength, but depends on wavelength; for this reason, the cross-correlation needs to be calculated in log$\lambda$ space.

Shifting images, and more generally, spatial sampling of images

The issue of shifting images is really an issue of spatial interpolation. The specific task of needing to shift images brings up a more fundamental discussion of spatial sampling of images; what knowledge is retained, and what is lost, when one observes only sampling of the spatial intensities?

There are two issues you want to consider when shifting images: you want to know if they will conserve flux and you want to know if they will preserve spatial information. Generally, we probably want to require flux conservation, and we want to preserve as much spatial information as necessary. You may also be interested in carrying along error information for each pixel values.

If you know the underlying shape, you can determine the value at every location. In reality, however, we sample or bin the underlying function. There is a distinction between sampling and binning. Binning becomes equivalent to sampling when your function in well approximated by a linear function across the pixel, but not if there is any curvature. In astronomical detectors, we always bin the underlying function. But we can also talk about sampling an underlying function that is smoothed by a pixel boxcar function, and hence about recovering the binned value at different pixel centerings.

For many imaging applications, there is a physical process (e.g., seeing) which limits the amount of high spatial frequency information. When it is the case that there is an upper frequency cutoff in the power spectrum of the objects being observed (band-limited), it is possible to recover the full function from samples of it if the samples are sufficiently fine. This is known as the sampling theorem. It says that if you have a band-limited function, you can recover the function if you have samples at spacing not exceeding 0.5/$\nu_{c}^{}$ , where $\nu_{c}^{}$ is the spatial frequency at which the power spectrum goes to zero. For example, if you say seeing wipes out scales less than 0.2 arcsec, you need to sample at 0.1 arcsec to recover the full function.

In practice, we don't measure a frequency cutoff and we don't rigorously define critical sampling, but we talk about it anyway. Generally, in astronomy, critical sampling is sampling at at least twice the full-width- half-maximum of the seeing disk (or the size of the Airy disk for diffraction limited applications). Note that this definition is pretty loose considering we have square, not circular pixels. Probably something like three samples is more reasonable.

If you can recover the underlying function, then you can sample it at different locations. You can do this using the sampling theorem by something known as sinc interpolation. This works by filtering by a box function in transform space, which is equivalent to convolution with a sinc function in real space.

However, sinc interpolation may not be accurate unless you are better than critically sampled. If you are undersampled, it fails miserably. Near critical sampling, it can lead to non-flux conserving interpolation, which you generally really want to avoid. So it is not preferred unless you are very well sampled.

Also, note that even for astronomical data that may be approximated as band-limited, there may be objects in the data that do not have the same spatial smearing, e.g. cosmic rays. These can turn into relatively nasty things after sinc interpolation!

Perhaps the simplest method of shifting is nearest pixel shifts, and in fact, it turns out to be a decent way to shift. The concept is simple: determine the shift between images, round it to the nearest integer, shift your image by that amount trivially, and add the results. This obviously conserves flux, and preserves full error information. The disadvantage is that it leads to some loss of spatial information (e.g. if you add an image with another that is off by 1/2 of a pixel), but this may not be so much of an issue, depending on the application.

More sophisticated routines are not necessarily better. Generally, these schemes use some sort of polynomial interpolation, where a nth order polynomial that runs through the n+1 points surrounding the desired location is determined. This is often known as Lagrangian interpolation. For n = 1 , this is just linear interpolation, or, in two dimensions, bilinear interpolation. In 2D, you can do bilinear interpolation one dimension at a time. The formula is just:

F(x, y₁) = F(x₁, y₁) + (x - x₁)*(F(x₂, y₁) - F(x₁, y₁)

F(x, y₂) = F(x₁, y₂) + (x - x₁)*(F(x₂, y₂) - F(x₁, y₂)

F(x, y) = F(x, y₁) + (y - y₁)*(F(x, y₂) - F(x, y₁)

This can be extended to higher order. Often, bicubic is used. However, higher order doesn't necessarily buy you anything. One has to be very careful with higher order schemes, because in their simplest implementations, they may not conserve flux. It is straightforward to write higher order routines which fit binned functions and thus preserve flux, but some implementations you find may not do this.

In the case of very well sampled images, all of these routines do very well. As you get less well sampled, however, you need to be more careful. In fact, nearest pixel shifting does as well or better as more sophisticated routines. Also, all other techniques lead to correlated errors between pixels.

Spline interpolation: guarantees continuity of derivatives. This rarely comes up in the handling of images, but often comes up in other interpolation applications.

If you have an instrument that produces undersampled images, it is possible to regain some additional spatial information by dithering exposures at subpixel spacings. Several algorithms have been developed for combining such sets of observations; a popular one is the DRIZZLE algorithm developed at STScI (see, e.g. Fruchter & Hook, PASP 114, 144, 2002, and web references at STScI). However, beware that this can be tricky to do properly (see, e.g., Lauer, PASP 111, 756, 1999) . If possible, it may be advisable to model sets of images with an underlying model constrained by the individual observations, rather than combining the observations.

Spectroscopic data reduction: MOVED TO LAB

1. normal CCD processing: overscan, bias, dark, preflash.

2. flat fielding. Note problem that dome flats have spectral energy distribution of light source. ``Flatten'' the flats in the wavelength direction to preserve error analysis, i.e. remove the large scale wavelength dependence, but preserve the pixel-to-pixel response variations. In the spatial direction, flat fielding is like imaging, but often the requirements on accuracy are less stringent.

3. wavelength calibration. Use arc lamps with known lines. Identify lines, determine line centers (centroid or fitting), and fit function to centers vs. wavelength.

4. flux calibration: correction for throughput as a function of wavelength. Not always required, e.g. if measuring strengths relative to nearby continuum. Spectrophotometric standards, e.g. Massey et al. ApJ 328, 315 (1988).

   If fluxing is performed, usually also want to correct for atmospheric extinction as a function of wavelength and airmass: use of mean extinction coefficients.

5. Object reduction: extracting object spectrum (``tracing'' the object) and sky spectrum. Aperture extraction vs. optimal extraction. Caveats: spectral curvature.

6. Advanced topics: nod and shuffle, atmospheric feature correction (esp in IR).