Cosmic Ray Detection Algorithm

Overview

Cosmic rays are high-energy particles that create spurious bright pixels in CCD images. Detection and flagging is critical for accurate spectroscopy. pyKOSMOS++ uses a scipy-based implementation of the L.A.Cosmic algorithm (van Dokkum 2001).

Method: L.A.Cosmic Algorithm

L.A.Cosmic detects cosmic rays by comparing image sharpness to a smoothed version:

Algorithm Steps:

  1. Create Laplacian image (edge-detection kernel):

    \[\begin{split}L = \\nabla^2 I = \\frac{\\partial^2 I}{\\partial x^2} + \\frac{\\partial^2 I}{\\partial y^2}\end{split}\]

  2. Estimate noise model:

    \[\begin{split}\\sigma = \\sqrt{\\sigma_{read}^2 + I / g}\end{split}\]

  3. Compute significance:

    \[\begin{split}S = \\frac{L}{\\sigma}\end{split}\]

  4. Apply contrast criterion:

    Pixel is cosmic ray if:

    • S > sigma_clip (typically 5σ)

    • AND pixel value > contrast × median of neighbors

  5. Iteratively clean and re-detect (max 5 iterations)

Parameters

sigma_clip: 5.0

Detection threshold in sigma above noise

contrast: 3.0

Minimum contrast with neighbors

max_iterations: 5

Maximum cleaning iterations

readnoise: 3.7 e⁻

CCD read noise (from config)

gain: 1.4 e⁻/ADU

CCD gain (from config)

Implementation

from scipy.ndimage import laplace, median_filter

def detect_cosmic_rays(
    data: np.ndarray,
    sigma_clip: float = 5.0,
    contrast: float = 3.0,
    max_iterations: int = 5
) -> np.ndarray:
    """
    Detect cosmic rays using L.A.Cosmic algorithm.

    Returns
    -------
    mask : ndarray (bool)
        True for cosmic ray pixels
    """
    mask = np.zeros(data.shape, dtype=bool)

    for iteration in range(max_iterations):
        # Laplacian edge detection
        laplacian = laplace(data)

        # Noise model
        noise = np.sqrt(readnoise**2 + data / gain)

        # Significance
        significance = laplacian / noise

        # Detect outliers
        med_neighbors = median_filter(data, size=3)
        contrast_check = data > contrast * med_neighbors

        new_cosmic = (significance > sigma_clip) & contrast_check

        if np.sum(new_cosmic) == 0:
            break  # Converged

        mask |= new_cosmic

        # Clean image for next iteration
        data[mask] = med_neighbors[mask]

    return mask

Performance

Typical Results:

  • Detection rate: 0.5-2% of pixels flagged

  • False positive rate: <0.1%

  • Processing time: 2-5 seconds per 2048×515 image

When to adjust parameters:

  • More aggressive detection (sigma_clip=4.0): More false positives, better for faint cosmic rays

  • Less aggressive (sigma_clip=6.0): Fewer false positives, may miss faint events

  • Contrast threshold (contrast=2.0-5.0): Lower = more detections

Integration with Pipeline

Cosmic ray detection occurs after flat correction, before trace detection:

  1. Load and calibrate science frame (bias subtract, flat correct)

  2. Detect cosmic rays and create mask

  3. Store mask in Spectrum2D object

  4. Apply mask during extraction (zero weight to flagged pixels)

Impact on Extraction

Flagged pixels receive zero weight during optimal extraction, preventing cosmic rays from contaminating the extracted 1D spectrum.

Without cosmic ray rejection:

  • Bright spikes in 1D spectrum

  • Artificial emission lines

  • Reduced SNR due to outliers

With cosmic ray rejection:

  • Clean 1D spectrum

  • Improved SNR (5-10% typical)

  • Accurate line flux measurements

Limitations

Cannot detect:

  • Cosmic rays aligned with spectral traces (< 1% of events)

  • Very faint cosmic rays (< 3σ above background)

  • Saturated cosmic rays (may create artifacts)

Workarounds:

  • Use multiple exposures with sigma-clipping combination

  • Inspect diagnostic plots for remaining artifacts

  • Manual masking if needed

References

See Also