Interactive use of cleanest functionality

This notebook describes how to handle the cleaning of cosmic-ray pixels when working with single or multiple equivalent exposures.

Warning: We are not trying here to demonstrate the usefulness of a particular detection algorithm nor the benefits of using a particular interpolation method to replace cosmic-ray pixels. We are simply showing the basic Python code required to use programmatically some of the functionality available in tea-cleanest.

Note: If you use this code, please cite Cardiel et al. 2026 and the references of the corresponding detection algorithms employed:

• L.A.Cosmic: see References at the bottom of this link
• PyCosmic: cite Husemann et al. 2012
• deepCR: cite Zhang et al. 2020
• Cosmic-CoNN: cite Xu et al. 2023

from astropy.io import fits
from datetime import datetime
import numpy as np
from scipy import ndimage

import teareduce as tea
from teareduce import cleanest

time_ini = datetime.now()

Import specific packages

The following code checks whether some specific packages devoted to detecting and cleaning cosmic rays have already been installed. This is not necessary if you have previously successfully installed and run tea-cleanest.

import importlib

# Map: pip package name -> Python import name
required = {
    "ccdproc": "ccdproc",
    "maskfill": "maskfill",     # change if the import name differs
    "deepCR": "deepCR",         # case sensitive; adjust if module is 'deepcr'
    "cosmic-conn": "cosmic_conn"  # hyphen in pip, underscore in import (example)
}

missing = []

for pip_name, import_name in required.items():
    print(f"Checking availability of package: {import_name}")
    try:
        importlib.import_module(import_name)
    except ModuleNotFoundError:
        missing.append((pip_name, import_name))

if missing:
    msg_lines = ["This notebook requires the following packages:"]
    for pip_name, import_name in missing:
        msg_lines.append(
            f" - import name: '{import_name}' (install with: pip install {pip_name})"
        )
    raise ModuleNotFoundError("\n".join(msg_lines))

# Import all packages into the global namespace
for import_name in required.values():
    print(f"Importing: {import_name}")
    globals()[import_name] = importlib.import_module(import_name)

Checking availability of package: ccdproc
Checking availability of package: maskfill
Checking availability of package: deepCR
Checking availability of package: cosmic_conn
Importing: ccdproc
Importing: maskfill
Importing: deepCR
Importing: cosmic_conn

try:
    import PyCosmic
except ModuleNotFoundError as e:
    raise ModuleNotFoundError(
        "The 'teareduce.cleanest' module requires the 'PyCosmic' package.\n"
        "Please install this module using:\n"
        "`pip install git+https://github.com/nicocardiel/PyCosmic.git@test`"        
    ) from e

print("Importing: PyCosmic")

Importing: PyCosmic

Example data

Download the example data files used in this notebook.

for fname in [
    'notebooks/cleanest/examplecr1.fits',
    'notebooks/cleanest/examplecr2.fits',
    'notebooks/cleanest/examplecr3.fits',
    'notebooks/cleanest/examplecr1_cleaned_with_surface.fits',
]:
    tea.get_cookbook_file(fname)

File examplecr1.fits already exists locally. Skipping download.
File examplecr2.fits already exists locally. Skipping download.
File examplecr3.fits already exists locally. Skipping download.
File examplecr1_cleaned_with_surface.fits already exists locally. Skipping download.

Detect CRs in a single image using L.A.Cosmic

The following lines of code show how to run the cleaning of the first example image used in Simple execution (single image). The results will be compared with a previously computed cleaned image, which serves as a test of the code.

# Read the input file
with fits.open('examplecr1.fits') as hdul:
    data = hdul[0].data

Running L.A.Cosmic only once

As explained in the documentation of tea-cleanest, the L.A.Cosmic algorithm can be run twice (run 1 and run 2) on the same image using different values for some parameters. Before dealing with this case, let’s illustrate how to run L.A.Cosmic simply once, using the same parameters used in run 1.

Here we are using the cleanest.lacosmicpad() function, which is a wrapper around the cosmicray_lacosmic() function from the ccdproc package. This function pads the input image array before applying L.A.Cosmic (this helps the algorithm locate CR pixels close to the image borders). If inbkg or invar arrays are provided, they are also padded accordingly. After processing, the padding is removed to return an array of the original size.

Users should consult the cosmicray_lacosmic() function documentation provided by the ccdproc package.

pad_width = 10

cleandata_lacosmic, mask_crfound = cleanest.lacosmicpad(
    pad_width=pad_width,
    ccd=data,
    gain=1.0,
    readnoise=6.5,
    sigclip=5.0,
    sigfrac=0.3,
    objlim=5.0,
    niter=4,
    fsmode='median',
    verbose=True,
)

Running L.A.Cosmic implementation from ccdproc version 2.4.3

(please wait...)

Starting 4 L.A.Cosmic iterations
Iteration 1:
6419 cosmic pixels this iteration
Iteration 2:
353 cosmic pixels this iteration
Iteration 3:
29 cosmic pixels this iteration
Iteration 4:
17 cosmic pixels this iteration

Done!

Check image shapes

data.shape, cleandata_lacosmic.shape, mask_crfound.shape

((1596, 2016), (1596, 2016), (1596, 2016))

Clean the detected CRs

Using the cosmic-ray pixel mask, we can proceed to interpolate the cosmic rays using the desired method.

The cosmic-ray pixels are initially dilated by the specified number of pixels (0 means no dilation). The resulting flagged pixels are grouped into cosmic-ray features. Each cosmic-ray feature is then interpolated using the specified interpolation method.

Note that the interpolation methods lacosmic and auxfile, available in the interactive use of tea-cleanest are not implemented in this function because both cases are simply an inmediate replacement of the cosmic ray pixels in the data array by the corresponding pixels in another array using the mask array. Therefore, these two methods do not require any interpolation algorithm.

In the following example we are replacing the detected pixels using the surface interpolation (interp_method='s').

data_cleaned, masked = cleanest.interpolate(data, mask_crfound, interp_method='s', npoints=2, degree=1, debug=True)

Number of cosmic ray pixels to be cleaned: 6590

Number of cosmic rays (grouped pixels)...:  749

Number of cosmic rays cleaned............:  749

For a more detailed explanation of the available possibilities, have a look at the documentation of this function.

?cleanest.interpolate

Signature:
cleanest.interpolate(
    data,
    mask_crfound,
    dilation=0,
    interp_method=None,
    npoints=None,
    degree=None,
    debug=False,
)
Docstring:
Interpolate pixels identified in a mask.

The original data and mask are not modified. A copy of both
arrays are created and returned with the interpolated pixels.

Cosmic-ray pixels are initially dilated by the specified number
of pixels. The resulting flagged pixels are grouped into cosmic-ray
features. Each cosmic-ray feature is then interpolated using the
specified interpolation method.

Note that the interpolation methods `lacosmic` and `auxfile`,
available in the interactive use of **tea-cleanest** are not
implemented in this function because both cases are simply
an inmediate replacement of the cosmic ray pixels in the data
array by the corresponding pixels in another array using the
mask array. Therefore, these two methods do not require any
interpolation algorithm.

Parameters
----------
data : 2D numpy.ndarray
    The image data array to be processed.
mask_crfound : 2D numpy.ndarray of bool
    A boolean mask array indicating which pixels are flagged
    and need to be interpolated (True = pixel to be fixed).
dilation : int, optional
    The number of pixels to dilate the masked pixels before
    interpolation.
interp_method : str, optional
    The interpolation method to use. Options are:
    'x' : Polynomial interpolation in the X direction using neighbor pixels
    'y' : Polynomial interpolation in the Y direction using neighbor pixels
    's' : Surface fit (degree 1) interpolation using neighbor pixels
    'd' : Median of neighbor pixels interpolation using neighbor pixels
    'm' : Mean of neighbor pixels interpolation using neighbor pixels
    'k' : Maskfill method, as described in van Dokkum & Pasha (2024)
npoints : int, optional
    The number of points to use for interpolation. This
    parameter is relevant for 'x', 'y', 's', 'd', and 'm' methods.
degree : int, optional
    The degree of the polynomial to fit. This parameter is
    relevant for 'x' and 'y' methods.
debug : bool, optional
    If True, print debug information and enable tqdm progress bar.

Returns
-------
cleaned_data : 2D numpy.ndarray
    The image data array with cosmic rays cleaned.
mask_fixed : 2D numpy.ndarray of bool
    The updated boolean mask array indicating which pixels
    have been fixed.

Notes
-----
This function has been created to clean cosmic rays without
the need of a GUI interaction. It can be used in scripts
or batch processing of images.
File:      ~/s/teareduce/z_devel/teareduce/src/teareduce/cleanest/interpolate.py
Type:      function

Running L.A.Cosmic twice

Although the prior execution of L.A.Cosmic has allowed the removal of many cosmic rays in the initial image, visual inspection reveals that the tails of some cosmic rays (i.e., pixels with weaker but still affected signals) have not been interpolated. The way tea-cleanest addresses this situation is by running L.A.Cosmic a second time, using a lower threshold to detect cosmic rays. To prevent the appearance of many false positives, the code groups suspected pixels into cosmic-ray features. The newly detected cosmic-ray features from the second run are only retained if they are in contact with cosmic-ray features already identified in the first execution of L.A.Cosmic.

cleandata_lacosmic2, mask_crfound2 = cleanest.lacosmicpad(
    pad_width=pad_width,
    ccd=data,
    gain=1.0,
    readnoise=6.5,
    sigclip=3.0,      # modified
    sigfrac=0.3,
    objlim=5.0,
    niter=4,
    fsmode='median',
    verbose=True,
)

Running L.A.Cosmic implementation from ccdproc version 2.4.3

(please wait...)

Starting 4 L.A.Cosmic iterations
Iteration 1:
9375 cosmic pixels this iteration
Iteration 2:
810 cosmic pixels this iteration
Iteration 3:
107 cosmic pixels this iteration
Iteration 4:
38 cosmic pixels this iteration

Done!

np.sum(mask_crfound), np.sum(mask_crfound2)

(np.int64(6590), np.int64(9954))

The merging of both masks, following the strategy previously mentioned, is carried out by the auxiliary function merge_peak_tail_masks().

mask_crfound = cleanest.merge_peak_tail_masks(
    mask_peaks=mask_crfound, 
    mask_tails=mask_crfound2,
    verbose=True
)

Number of cosmic ray pixels (peaks)...............: 6590
Number of cosmic ray pixels (tails)...............: 9954
Number of cosmic ray pixels (merged peaks & tails): 7648

We then proceed to interpolate the cosmic-ray pixels using the new mask.

data_cleaned, masked = cleanest.interpolate(
    data=data, 
    mask_crfound=mask_crfound, 
    interp_method='s', 
    npoints=2, 
    degree=1, 
    debug=True
)

Number of cosmic ray pixels to be cleaned: 7648

Number of cosmic rays (grouped pixels)...:  748

Number of cosmic rays cleaned............:  748

Comparison with reference image

One of the downloaded images at the beginning of the execution of this notebook is the file examplecr1_cleaned_with_surface.fits, which was generated interactively using tea-cleanest. We can compare this file with the result from the last section.

reference_filename = "examplecr1_cleaned_with_surface.fits"
with fits.open(reference_filename) as hdul:
    reference_data = hdul[0].data
    reference_mask = hdul["CRMASK"].data

Check data array

np.array_equal(data_cleaned, reference_data)

True

Check CR mask array

np.array_equal(masked, reference_mask)

True

Using other CR detection algorithms for single images

You can choose any suitable detection algorithm before making use of the cleanest.interpolation() function. Let’s see how to do it using other algorithms incorporated in tea-cleanest.

PyCosmic

In this case one can proceed by calling directly the det_cosmics() function available in the PyCosmic package.

?PyCosmic.det_cosmics

Signature:
PyCosmic.det_cosmics(
    data,
    sigma_det=5,
    rlim=1.2,
    iterations=5,
    fwhm_gauss=[2.0, 2.0],
    replace_box=[5, 5],
    replace_error=1000000.0,
    increase_radius=0,
    gain=1.0,
    rdnoise=1.0,
    bias=0.0,
    verbose=False,
)
Docstring:
Detects and removes cosmics from astronomical images based on Laplacian edge
detection scheme combined with a PSF convolution approach.

IMPORTANT:
The image and the readout noise are assumed to be in units of electrons.
The image also needs to be BIAS subtracted! The gain can be entered to convert the image from ADUs to
electros, when this is down already set gain=1.0 as the default. If ncessary a homegnous bias level can
be subtracted if necessary but default is 0.0.

 Parameters
 --------------
 data: ndarray
         Two-dimensional array representing the input image in which cosmic rays are detected.
 sigma_det: float, default: 5.0
         Detection limit of edge pixel above the noise in (sigma units) to be detected as comiscs
 rlim: float, default: 1.2
         Detection threshold between Laplacian edged and Gaussian smoothed image
 iterations: integer, default: 5
         Number of iterations. Should be >1 to fully detect extended cosmics
 fwhm_gauss: list of floats, default: [2.0, 2.0]
         FWHM of the Gaussian smoothing kernel in x and y direction on the CCD
 replace_box: list integers, default: [5,5]
         median box size in x and y to estimate replacement values from valid pixels
 replace_error: float, default: 1e6
         Error value for bad pixels in the comupted error image
 increase_radius: integer, default: 0
         Increase the boundary of each detected cosmic ray pixel by the given number of pixels.
 gain: float, default=1.0
         Value of the gain in units of electrons/ADUs
 rdnoise: float, default=1.0
         Value of the readout noise in electrons
 bias: float, default=0.0
         Optional subtraction of a bias level.
 verbose: boolean, default: False
         Flag for providing information during the processing on the command line

 Ouput
 -------------
 out: Image class instance
     Result of the detection process is an Image which contains .data, .error, .mask as attributes for the
     cleaned image, the internally computed error image and a mask image with flags for cosmic ray pixels.

 Reference
 --------------
 Husemann et al. 2012, A&A, Volume 545, A137 (https://ui.adsabs.harvard.edu/abs/2012A%26A...545A.137)
File:      ~/venv_tea/lib/python3.12/site-packages/PyCosmic/det_cosmics.py
Type:      function

out = PyCosmic.det_cosmics(
    data=data,
    sigma_det=5.0,
    rlim=1.2,
    iterations=5,
    fwhm_gauss=[2.5, 2.5],
    replace_box=[5, 5],
    replace_error=1e6,
    increase_radius=0,
    gain=1.0,
    rdnoise=6.5,
    bias=0,
    verbose=True,
)

A value of 6.500000 is used for the electron read-out noise.
Start the detection process using.
Start iteration 1
Total number of detected cosmics: 2966 out of 3217536 pixels
Start iteration 2
Total number of detected cosmics: 5617 out of 3217536 pixels
Start iteration 3
Total number of detected cosmics: 6333 out of 3217536 pixels
Start iteration 4
Total number of detected cosmics: 6414 out of 3217536 pixels
Start iteration 5
Total number of detected cosmics: 6418 out of 3217536 pixels

Although the previous function already returns a cleaned version of the image in out.data, the use of cleanest.interpolation() allows you to use simply the mask stored in out.mask and select any of the available interpolation algorithms in the latter function.

mask_crfound = out.mask.astype(bool)

data_cleaned, masked = cleanest.interpolate(
    data=data, 
    mask_crfound=mask_crfound,
    interp_method='s', 
    npoints=2, 
    degree=1, 
    debug=True
)

Number of cosmic ray pixels to be cleaned: 6418

Number of cosmic rays (grouped pixels)...:  806

Number of cosmic rays cleaned............:  806

deepCR

Here we make use of two functions from the deepCR package

?deepCR.deepCR

Init signature: deepCR.deepCR(mask='ACS-WFC', inpaint=None, device='CPU', hidden=32)
Docstring:      <no docstring>
Init docstring:
    Instantiation of deepCR with specified model configurations

Parameters
----------
mask : str
    Either name of existing deepCR-mask model, or file path of your own model (incl. '.pth')
inpaint : (optional) str
    Name of existing inpainting model to use. If left as None then by default use a simple 5x5 median mask
    sampling for inpainting
device : str
    One of 'CPU' or 'GPU'
hidden : int
    Number of hidden channel for first deepCR-mask layer. Specify only if using custom deepCR-mask model.
Returns
-------
None
    
File:           ~/venv_tea/lib/python3.12/site-packages/deepCR/model.py
Type:           type
Subclasses:     

# Initialize the deepCR model
mdl = deepCR.deepCR(mask="ACS-WFC")

?mdl.clean

Signature:
mdl.clean(
    img0,
    threshold=0.5,
    inpaint=False,
    binary=True,
    segment=True,
    patch=1024,
    n_jobs=1,
)
Docstring:
    Identify cosmic rays in an input image, and (optionally) inpaint with the predicted cosmic ray mask
:param img0: (np.ndarray) 2D input image conforming to model requirements. For HST ACS/WFC, must be from
_flc.fits and in units of electrons in native resolution.
:param threshold: (float; [0, 1]) applied to probabilistic mask to generate binary mask
:param inpaint: (bool) return clean, inpainted image only if True
:param binary: return binary CR mask if True. probabilistic mask if False
:param segment: (bool) if True, segment input image into chunks of patch * patch before performing CR rejection.
  Used for memory control.
:param patch: (int) Use 256 unless otherwise required. if segment==True, segment image into chunks of
  patch * patch.
:param n_jobs: (int) number of jobs to run in parallel, passed to `joblib.` default: 1.
:return: CR mask and (optionally) clean inpainted image
File:      ~/venv_tea/lib/python3.12/site-packages/deepCR/model.py
Type:      method

mask_crfound, cleandata_deepcr = mdl.clean(
    img0=data,
    threshold=0.5,
    inpaint=True,
)

Although the previous function already returns a cleaned version of the image in cleandata_deepcr, the use of cleanest.interpolation() allows you to use simply the mask stored in mask_crfound and select any of the available interpolation algorithms in the latter function.

data_cleaned, masked = cleanest.interpolate(
    data=data, 
    mask_crfound=mask_crfound,
    interp_method='s', 
    npoints=2, 
    degree=1, 
    debug=True
)

Number of cosmic ray pixels to be cleaned: 7906.0

Number of cosmic rays (grouped pixels)...:    731

Number of cosmic rays cleaned............:    731

Cosmic-CoNN

Here we use the functionality available in the cosmic_conn package.

# Initialize the generic ground-imaging model
cr_model = cosmic_conn.init_model("ground_imaging")

/Users/cardiel/venv_tea/lib/python3.12/site-packages/torch/__init__.py:1275: UserWarning: torch.set_default_tensor_type() is deprecated as of PyTorch 2.1, please use torch.set_default_dtype() and torch.set_default_device() as alternatives. (Triggered internally at /Users/runner/work/pytorch/pytorch/pytorch/torch/csrc/tensor/python_tensor.cpp:436.)
  _C._set_default_tensor_type(t)

If we get a warning message here is because the installed cosmic_conn package is using an old PyTorch function (torch.set_default_tensor_type()) that has been deprecated in PyTorch 2.1+. The PyTorch developers want people to use the newer functions torch.set_default_dtype() and torch.set_default_device() instead.

# The model outputs a CR probability map
print(f"Running Cosmic-CoNN version: {cosmic_conn.__version__}  (please wait...)")
cr_prob = cr_model.detect_cr(data.astype(np.float32))
print("Done!")

Running Cosmic-CoNN version: 0.5.0  (please wait...)
Done!

# Set probability threshold
threshold = 0.5

mask_crfound = cr_prob > threshold

In this case we only have a detection mask and not a cleaned version of the data. To obtain a cleaned version, we need to use cleanest.interpolate().

data_cleaned, masked = cleanest.interpolate(
    data=data, 
    mask_crfound=mask_crfound,
    interp_method='s', 
    npoints=2, 
    degree=1, 
    debug=True
)

Number of cosmic ray pixels to be cleaned: 6858

Number of cosmic rays (grouped pixels)...:  830

Number of cosmic rays cleaned............:  830

Detect CRs in multiple equivalent exposures

When several equivalent exposures are available, one can try to clean each of them making use of the information from the remaining images. The simplest situation is considering only two exposures.

When more than two images are available, more possibilities can be considered. For such cases, the auxiliary function cleanest.combine_arrays() can be useful (see below).

Simple case: two equivalent exposures

When considering only two exposures, one can detect CRs in one image and replace the suspected pixels using information from the second image. The roles of both images can be interchanged to clean the second image using information from the first exposure. The basic assumption behind this strategy is that the same pixels have not been hit by cosmic rays in both exposures, which may not be the case for long exposure times.

Let’s illustrate this approach cleaning the image in examplecr1.fits using the data in examplecr2.fits.

with fits.open('examplecr1.fits') as hdul:
    data1 = hdul[0].data

with fits.open('examplecr2.fits') as hdul:
    data2 = hdul[0].data

We are using the L.A.Cosmic procedure, which is relatively fast.

The underscore _ in the next cell is used to discard the first return value (the cleaned array) from lacosmicpad(), following the Python convention for intentionally unused variables. Only the cosmic ray mask is needed for subsequent processing.

list_sigclip = [5.0, 3.0]   # run L.A.Cosmic twice
list_mask_crfound = []

for sigclip in list_sigclip:
    print(f"Using L.A.Cosmic with {sigclip=}")
    _, mask_crfound = cleanest.lacosmicpad(
        pad_width=10,
        ccd=data1,
        gain=1.0,
        readnoise=6.5,
        sigclip=sigclip,
        sigfrac=0.3,
        objlim=5.0,
        niter=4,
        fsmode='median',
        verbose=False,
    )
    list_mask_crfound.append(mask_crfound)

Using L.A.Cosmic with sigclip=5.0
Using L.A.Cosmic with sigclip=3.0

mask_crfound = cleanest.merge_peak_tail_masks(
    mask_peaks=list_mask_crfound[0],
    mask_tails=list_mask_crfound[1],
    verbose=True
)

Number of cosmic ray pixels (peaks)...............: 6590
Number of cosmic ray pixels (tails)...............: 9954
Number of cosmic ray pixels (merged peaks & tails): 7648

Replace masked pixels in data1 using the information from data2.

data1_cleaned = data1.copy()
data1_cleaned[mask_crfound] = data2[mask_crfound]

As previously mentioned, once the first image has been cleaned, one can swap the roles of the two images to clean the second one making use of the first exposure.

To minimize the need to code all the previous steps, you can consider using the function cleanest.combine_arrays() described next.

General case

When the number of individual equivalent exposures is larger than 2, several possibilities can be considered to generate a combined cleaned version of all of them. To help in this task, the cleanest module offers a handy function combine_arrays(). Note that this function can also be used even in the case of having only 2 equivalent exposures.

Note: This function can also be used even in the case of having only 2 equivalent exposures.

The function combine_arrays() takes a list of 2D arrays, applies the specified cleaning method to each individual array using the chosen strategy, and then combines the cleaned arrays employing one of the available combination methods.

?cleanest.combine_arrays

Signature:
cleanest.combine_arrays(
    list_arrays,
    detection_algorithm,
    cleaning_strategy,
    combination_method,
    show_progress=False,
    return_array_mask_lists=False,
    **kwargs,
)
Docstring:
Combine arrays with different interpolation methods.

The function takes a list of 2D arrays, applies the specified
cleaning method to each individual array using the
chosen strategy, and then combines the cleaned arrays employing
one of the available combination methods.

Parameters
----------
list_arrays : list of 2D numpy.ndarray
    A list of 2D arrays to be combined.
detection_algorithm : str
    The algorithm used to detect cosmic rays. Supported algorithms are:
    - "lacosmic": Use the L.A.Cosmic algorithm for CR detection.
    - "pycosmic": Use the PyCosmic algorithm for CR detection.
    - "deepcr": Use the DeepCR algorithm for CR detection.
    - "conn": Use the Cosmic-CoNN algorithm for CR detection.
cleaning_strategy : str
    The strategy used for cleaning the arrays. Supported strategies are:
    - "crmethod": replace the cosmic ray pixels in each array with the values
      obtained by applying the specified detection algorithm to that array
      alone. Note that this strategy does not use the information from the
      other arrays in `list_arrays` to clean the cosmic ray pixels,
      but only to combine the results after cleaning.
    - "median_all": replace the cosmic ray pixels with the median of all
      the arrays in `list_arrays`.
    - "median_others": replace the cosmic ray pixels with the median of all the
        arrays in `list_arrays` except the one being cleaned.
    - "mean_others": replace the cosmic ray pixels with the mean of all the
        arrays in `list_arrays` except the one being cleaned.
    - "none": do not attempt to clean the arrays, just generate
        the masks of the detected cosmic rays and combine them using the
        masked version of the specified combination method.
combination_method : str
    The method used to combine the arrays. Supported methods are:
    - "median": median of the individually cleaned or masked arrays.
    - "mean": mean of the individually cleaned or masked arrays.
show_progress : bool, optional
    If True, print show_progress information during the combination process.
return_array_mask_lists : bool, optional
    If True, return a list of the individually cleaned arrays and another
    list with the corresponding masks.
**kwargs : dict
    Additional keyword arguments passed to the CR detection algorithm.

Returns
-------
combined_array : 2D numpy.ndarray
    The resulting array after combining the input arrays using
    the specified methods.
out_list_arrays: list of 2D numpy.ndarray, optional
    If return_list_arrays_masks is True, a list of the individually
    cleaned arrays.
out_list_masks: list of 2D numpy.ndarray of bool, optional
    If return_list_arrays_masks is True, a list of the corresponding
    masks of detected CRs for each array.
File:      ~/s/teareduce/z_devel/teareduce/src/teareduce/cleanest/combine_arrays.py
Type:      function

Note that the function first looks for cosmic rays in each individual exposure, generating a cleaned version and the corresponding mask for each of them. For that purpose, any of the 4 detection algorithms implemented in tea-cleanest can be used: lacosmic, pycosmic, deepcr, and conn. The corresponding algorithm is specified in the parameter detection_algorithm.

The user should define the cleaning_strategy:

• crmethod: replace the cosmic-ray pixels in each array with the values obtained by applying the specified detection algorithm to that array alone. Note that this strategy does not use the information from the other arrays in list_arrays to clean the cosmic-ray pixels, but only to combine the results after cleaning. Note that this option is not available when detection_algorithm=conn.

• median_all: replace the cosmic-ray pixels with the median of all the arrays in list_arrays in the corresponding pixels.

• median_others: replace the cosmic-ray pixels with the median of all the arrays in list_arrays except the one being cleaned.

• mean_others: replace the cosmic-ray pixels with the mean of all the arrays in list_arrays except the one being cleaned.

• none: do not attempt to clean the arrays, just generate the masks of the detected cosmic rays and combine them using the masked version of the specified combination method.

A combination method for the cleaned individual images should also be specified in combination_method. The two possible values are mean and median.

The show_progress boolean parameter indicates whether the function should provide a short description of the processing being carried out.

By setting the boolean optional parameter return_array_mask_lists to True, the function returns two additional lists: one containing the individually cleaned arrays and another one providing the corresponding pixel masks.

Finally, the user should include in this function all the parameters required by the chosen CR detection algorithm, excluding the parameter indicating the array to be cleaned, which will be defined properly by the combine_arrays() function as required when looping through the list_arrays. In particular, the parameter that must not be included is:

• ccd: for detection_algorithm=lacosmic

• data: for detection_algorithm=pycosmic

• img0: for detection_algorithm=deepcr

• image: for detection_algorithm=conn

Note: The function combine_arrays() can also execute the selected CR detection algorithm twice. This is the case if any of the input parameters for the corresponding method is a list of two values instead of a single one. This approach is used to detect the brightest CR pixels and the fainter CR tails with different thresholds. The resulting masks are properly merged in the final CR mask.

Using L.A.Cosmic

Let’s have a look at combine_arrays() in action, using L.A.Cosmic first.

with fits.open('examplecr1.fits') as hdul:
    data1 = hdul[0].data

with fits.open('examplecr2.fits') as hdul:
    data2 = hdul[0].data

with fits.open('examplecr3.fits') as hdul:
    data3 = hdul[0].data

Note that by setting sigclip=[5.0, 3.0], the detection algorithm will be run twice.

combined_array = cleanest.combine_arrays(
    list_arrays=[data1, data2, data3],
    detection_algorithm="lacosmic",  # lacosmic | pycosmic | deepcr | conn
    cleaning_strategy="crmethod",    # crmethod | median_all | median_other | mean_other | none
    combination_method="mean",       # median | mean
    show_progress=True,
    # parameters for the CR detection algorithm
    pad_width=10,
    gain=1.0,
    readnoise=6.5,
    sigclip=[5.0, 3.0],
    sigfrac=0.3,
    objlim=5.0,
    niter=4,
    fsmode='median',
    verbose=False,
)

Processing array 1/3...
Running L.A.Cosmic algorithm (run 1/2)...
Running L.A.Cosmic algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............: 6590
Number of cosmic ray pixels (tails)...............: 9954
Number of cosmic ray pixels (merged peaks & tails): 7648

Processing array 2/3...
Running L.A.Cosmic algorithm (run 1/2)...
Running L.A.Cosmic algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............:  7661
Number of cosmic ray pixels (tails)...............: 11106
Number of cosmic ray pixels (merged peaks & tails):  8833

Processing array 3/3...
Running L.A.Cosmic algorithm (run 1/2)...
Running L.A.Cosmic algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............:  6859
Number of cosmic ray pixels (tails)...............: 10127
Number of cosmic ray pixels (merged peaks & tails):  7925
Combining arrays using method 'mean'...

Combination complete!

The procedure for the other methods is similar.

Using PyCosmic

Note that by setting sigma_det=[5.0, 3.0], the detection algorithm will be run twice.

In this sense, it is important to note that the fwhm_gauss and replace_box parameters are already defined as a list of two values, so if you need to use different values in the two possible executions of the algorithm, these parameters must be set as a list of two-element lists.

combined_array = cleanest.combine_arrays(
    list_arrays=[data1, data2, data3],
    detection_algorithm="pycosmic",  # lacosmic | pycosmic | deepcr | conn
    cleaning_strategy="crmethod",    # crmethod | median_all | median_other | mean_other | none
    combination_method="mean",       # median | mean
    show_progress=True,
    # parameters for the CR detection algorithm
    sigma_det=[5.0, 3.0],
    rlim=1.2,
    iterations=5,
    fwhm_gauss=[2.5, 2.5],
    replace_box=[5, 5],
    replace_error=1e6,
    increase_radius=0,
    gain=1.0,
    rdnoise=6.5,
    bias=0,
    verbose=False,
)

Processing array 1/3...
Running PyCosmic algorithm (run 1/2)...
Running PyCosmic algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............: 6418
Number of cosmic ray pixels (tails)...............: 8958
Number of cosmic ray pixels (merged peaks & tails): 7242

Processing array 2/3...
Running PyCosmic algorithm (run 1/2)...
Running PyCosmic algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............:  7464
Number of cosmic ray pixels (tails)...............: 10031
Number of cosmic ray pixels (merged peaks & tails):  8394

Processing array 3/3...
Running PyCosmic algorithm (run 1/2)...
Running PyCosmic algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............: 6668
Number of cosmic ray pixels (tails)...............: 9185
Number of cosmic ray pixels (merged peaks & tails): 7516
Combining arrays using method 'mean'...

Combination complete!

Using deepCR

combined_array = cleanest.combine_arrays(
    list_arrays=[data1, data2, data3],
    detection_algorithm="deepcr",    # lacosmic | pycosmic | deepcr | conn
    cleaning_strategy="crmethod",    # crmethod | median_all | median_other | mean_other | none
    combination_method="mean",       # median | mean
    show_progress=True,
    # parameters for the CR detection algorithm
    threshold=[0.7, 0.3],
    inpaint=True,
)

Processing array 1/3...
Running DeepCR algorithm (run 1/2)...
Running DeepCR algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............: 7457
Number of cosmic ray pixels (tails)...............: 8355
Number of cosmic ray pixels (merged peaks & tails): 8347

Processing array 2/3...
Running DeepCR algorithm (run 1/2)...
Running DeepCR algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............: 8663
Number of cosmic ray pixels (tails)...............: 9686
Number of cosmic ray pixels (merged peaks & tails): 9679

Processing array 3/3...
Running DeepCR algorithm (run 1/2)...
Running DeepCR algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............: 7658
Number of cosmic ray pixels (tails)...............: 8606
Number of cosmic ray pixels (merged peaks & tails): 8595
Combining arrays using method 'mean'...

Combination complete!

Using Cosmic-CoNN

In this case, the cleaning_strategy cannot be crmethod since Cosmic-CoNN does not provide a cleaned array, but only a mask of detected CRs.

combined_array = cleanest.combine_arrays(
    list_arrays=[data1, data2, data3],
    detection_algorithm="conn",      # lacosmic | pycosmic | deepcr | conn
    cleaning_strategy="median_others",  # crmethod | median_all | median_other | mean_other | none
    combination_method="mean",       # median | mean
    show_progress=True,
    # parameters for the CR detection algorithm
    threshold=[0.7, 0.3],
)

Processing array 1/3...
Running Cosmic-CoNN algorithm (run 1/2)...
Running Cosmic-CoNN algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............: 6231
Number of cosmic ray pixels (tails)...............: 7674
Number of cosmic ray pixels (merged peaks & tails): 7152

Processing array 2/3...
Running Cosmic-CoNN algorithm (run 1/2)...
Running Cosmic-CoNN algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............: 7363
Number of cosmic ray pixels (tails)...............: 8670
Number of cosmic ray pixels (merged peaks & tails): 8328

Processing array 3/3...
Running Cosmic-CoNN algorithm (run 1/2)...
Running Cosmic-CoNN algorithm (run 2/2)...
Merging masks of detected cosmic rays from both runs...
Number of cosmic ray pixels (peaks)...............: 6667
Number of cosmic ray pixels (tails)...............: 8191
Number of cosmic ray pixels (merged peaks & tails): 7628
Combining arrays using method 'mean'...

Combination complete!

tea.elapsed_time_since(time_ini)

system...........: Darwin
release..........: 25.2.0
machine..........: arm64
node.............: Nicolass-MacBook-Pro.local
Python executable: /Users/cardiel/venv_tea/bin/python3.12
teareduce version: 0.7.1
Initial time.....: 2026-02-09 16:01:50.600552
Final time.......: 2026-02-09 16:07:39.169698
Elapsed time.....: 0:05:48.569146