IDL code for (1)-(2) is contained in rowdemo.pro.
IDL code for (3)-(5) is contained in gain.pro.
This analysis represents the preprocessing of Venus image data taken with the SpeX guider at the IRTF on July 6, 2004. The raw data for all examples below is contained in the file images/data1136.a.fits. After input, the raw data is first clipped to 0 and above. A simple mask is then created from all the positive pixel values. It appears that all negative pixel values in the file represent bad pixels, mostly due to CCD saturation and wraparound.
Row median values are computed from the data and mask with the function rowmed. Masked pixels are not included in the median calculation. It can be seen that the row medians oscillate around the curve for their true values, and that their Fourier transform contains a peak at the maximum spatial frequency. This oscillation is called picket noise, since it resembles the true median curve modulated by a square wave. The visual effect of this noise is a repeated brightening and dimming of alternating rows of pixels in the data.
The row median spectrum is lowpass filtered to about the second zero (arbitrarily chosen) and inverted to produce an approximation of the true median curve. Note that the effect of Fourier filtering introduces edge effects (similar to Gibbs phenomenon) at both ends of the approximation.
A row gain is calculated as the row medians divided by the filtered median curve. This gain oscillates about 1.0, and represents the departure of the observed medians from the actual medians. The data is then demodulated row by row using the function fixrows, which divides each row by the gain. Note that if the medians had been highpass filtered instead of lowpass filtered, this gain would oscillate about 0.0, and demodulation could not be performed in this way.
The demodulated data then has row medians equal to the desired curve (including the Fourier edge effects).
Closeups of the raw and filtered data show that the picket pattern has been removed from the background and the clouds, but not from the edge of the saturated limb. Nevertheless, after the picket noise has been removed, a faint embossed pattern which is more complicated than the original noise remains in the background.
To create a mask which rejects the saturated pixels in the limb of the planet, first the simple positive-value mask from (1) is median filtered with a width of 5.
In order to remove the additional picket noise at the edge of the limb, a horizontal edge-detection filter of size 5x3 pixels is created. The filter is normalized so that the sum of each column is equal to 0.
IDL> print, filt -0.500000 -0.500000 -0.500000 -0.500000 -0.500000 1.00000 1.00000 1.00000 1.00000 1.00000 -0.500000 -0.500000 -0.500000 -0.500000 -0.500000
The row-median-filtered data is then convolved by this edge filter to produce a response that hilites the edges of the saturated meniscus.
The filter response is used to create an additional mask in which good pixels are within less than 2.5 standard deviations of the mean of the response.
The two masks are then combined.
A brute-force 25-iteration median-filtering of width 5 is used to fill in the holes in the combined mask. The filtered data is then multiplied by the combined mask to remove the saturation. All remaining pixel values are now within the valid output range of the CCD, including those values at the edge of the limb. The mask can be extended at its edges to remove additional pixels, but it does not fit exactly or smoothly into the cusps of the limb.
To calculate an estimate of the CCD ADU-to-exposure relationship for each pixel, 36 frames of sky data were taken at 3 different timed exposures. These files include frames 955-990 of the July 6 data. The vector of exposure times (in seconds) for all 36 frames are
Exposure = 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 10.0000 20.0000 20.0000 20.0000 20.0000 20.0000
At each pixel, 36 (x, y) points were fitted to a line with the IDL function linfit, where the x values contain the vector of exposure times for each sky frame. 3 output arrays were computed for the parameters a and b in
and for the chi^2 values of each fit.
The values as initially computed are dominated by extremes for bad pixels.
max(a), min(a), mean(a), median(a), stddev(a) 315.486 -963.174 -42.1679 -41.9709 17.2903 max(b), min(b), mean(b), median(b), stddev(b) 13768.6 -5.33383 1002.92 1034.91 195.451 max(c2), min(c2), mean(c2), median(c2), stddev(c2) 4.85170e+07 0.00000 174253. 146895. 237339.
To create a pixel mask which rejects the extremes, only those values of a, b, and chi^2 which were within 2.5 standard deviations of their mean were retained. The values of a, b, and chi^2 are then multipled by this mask.
max(a), min(a), mean(a), median(a), stddev(a) 1.05208 -85.3889 -40.7687 -41.5291 17.2945 max(b), min(b), mean(b), median(b), stddev(b) 1484.34 0.00000 983.510 1032.94 238.436 max(c2), min(c2), mean(c2), median(c2), stddev(c2) 767371. 0.00000 162903. 144169. 103610.
Note the following:
The following plots show examples of the best-fit lines for the pixels (200, 100:109), which are within the domain of the 2.5 std mask.
Most of the Venus planetary image frames were exposed at 1 second (4 coadds of 0.25 sec each). A flat for this 1 sec exposure was calculated as:
which is the estimated density of a 1 sec exposure at each linearly fitted pixel. This flat was then divided by its median value.
max(flat), min(flat), mean(flat), median(flat), stddev(flat) 1.44295 0.00000 0.951930 1.00000 0.230552
To remove the picket noise without using a Fourier method, the raw data was divided by the flat image at each good point in the 2.5 std mask. This effectively removes the alternating variations in row gain, but not the saturation noise at the edge of the limb.
Images stretched to show 10 cycles of wrapped values from min to max further demonstrate the elimination of the picket noise from the background and clouds. Flat-fielding does not, however, romove all of the spatially correlated abberations in the CCD array, as can be seen by the vertical, and other, 'ripples' in the data. This may be because a dark frame (with negative intercept?) was not subtracted from the data before division by the flat frame.
The row medians as calculated by rowmed show that the observed values now match the expected values as in (1), but without the Fourier edge effects. However, there is still some oscillation about the true value due to the saturation noise along the limb, especially for row numbers less than 100.
Masking the flat-fielded data by the satuartion mask calculated in (2), and then calculating the row medians produces more accurate values for those rows which intersect the limb. Both the oscillations, and the ambient values, of these rows are reduced.
Additional sky frames, made at many different exposure times between zero and the saturation time of the CCD array, could be used to fit a higher order curve to the ADU-to-exposure relationship for each pixel. This would in turn allow a more accurate calculation of a flat field (and dark frame?) than is possible with a straight line, and might remove more of the spatially correlated abberations.