Finding and Removing Cosmic Rays

Finding Cosmic Rays

We look for cosmic rays using a sequence of tests, each more discriminating (and more time-consuming) than the preceeding. Initially the frame is searched for pixels that appear to be part of a cosmic ray event; then all of these contaminated pixels are assembled into cosmic ray events, and the pixels around these are examined for evidence of contamination, using a lower standard of proof.

In the following descriptions, C is the candidate cosmic ray pixel, and other pixels are referred to as (e.g.) NW.

               -1     0    1
             1 NW     N    NE
             0  W     C    E
            -1 SW     S    SE

Searching for contaminated pixels follows the following sequence:

If the code is compiled suitably, pixels that are labelled as interpolated, saturated, or not checked are ignored. This is usually not done, as it makes it hard to find CRs next to bad columns and bleed trails (Exercise for the reader: consider two equally bright CR pixels, one above the other, adjacent to a bad column. Why is the CR finder unable to find them?).
Then compare C with the four mean values of pairs of neighbouring pixels, (W--E), (N--S), (SW--NE), (NW--SE); if C exceeds any one of these by more than the min_sigma*sigma_sky the pixel is taken to be a cosmic ray candidate.
Next, see if the data is within the band limit, i.e. we know that no astronomical source can be sharper than the PSF, so the ratio of a pixel's value to that of a pixel a distance d away cannot exceed the ratio of the central value of a star to the value of a pixel at the same distance, d.
Allowing for noise, this condition becomes
```
(C - F*N(C) - sky) > F1*PSF(d)*(mean - F*N(mean) - sky)
```
where PSF(d) is the value of the PSF at a distance d, mean is the average of two pixels each a distance d away, and N(p) is p's standard deviation. F is an input parameter (called cond3_fac), with default value 3.0; F1 is a fiddle factor for the noise (it's called cond3_fac2 and has default value 0.6), used to deter photo from labelling the centres of bright stars as cosmic rays.
This inequality is evaluated for the four pairs N--S (d==1), E--W (d==1), SW--NE (d = sqrt(2)), and NW--SE (d = sqrt(2)), and if any satisfy it, we've got a cosmic ray.

After a contaminated pixel has been detected, its value is temporarily replaced by the lowest of the means determined in the last step listed above, in order to make it easier to find further pixels contaminated by the same cosmic ray.

When the whole frame has been processed, contiguous contaminated pixels are merged into single cosmic rays. If the peak pixel is below some minimum (cr_min_e) number of electrons above the sky level, the candidate cosmic ray is rejected, otherwise it's added to the MASK_CR mask.

Next, all cosmic ray pixels are removed.

All the pixels adjacent to each cosmic ray are then examined, using the above set of tests, but with the min_sigma and cond3_fac values relaxed by 50%. Any contaminated pixels are added to their `parent' cosmic ray.

Finally, if the keep flag is set, the original values of all cosmic ray pixels are reinstated; otherwise the MASK_INTERP bit is set for all interpolated pixels

Removing Cosmic Rays

Once all contaminated pixels have been found (as described above), they are interpolated over using a four-point interpolation scheme as described elsewhere (e.g. horizontally we calculate

	-0.2737*I[-2] + 0.7737*I[-1] + 0.7737*I[1] - 0.2737*I[2]

as a possible replacement for I[0]). We interpolate in all four directions (e.g. SW--NE, W--E) and choose the minimum value, after rejecting any involving bad (MASK_INTERP or MASK_NOTCHECKED) pixels. If none of the interpolated values is acceptable, the pixel is replaced by the sky value.

Because we have chosen the minimum of a number of estimates, the resulting value is biased a little low. We correct for this by assuming that there were only two uncontaminated values, and thus add 0.56*(sky sigma).

Restrictions on CR code

Cosmic ray contaminated pixels in the first and last row or column are more likely to be missed.
No more than 10000 cosmic ray contaminated pixels can be found in any one frame.