pypeit.core.procimg module

Module for image processing core methods

pypeit.core.procimg.base_variance(rn_var, darkcurr=None, exptime=None, proc_var=None, count_scale=None)[source]

Calculate the “base-level” variance in a processed image driven by the detector properties and the additive noise from the image processing steps.

The full variance model (see variance_model()), \(V\), is:

\[V = s^2\ \left[ {\rm max}(0, C) + D t_{\rm exp} / 3600 + V_{\rm rn} + V_{\rm proc} \right] + \epsilon^2 {\rm max}(0, c)^2\]

where:

\(c=s\ C\) are the rescaled observed sky + object counts,

\(C\) is the observed number of sky + object counts,

\(s=s\prime / N_{\rm frames}\) is a scale factor derived from the (inverse of the) flat-field frames plus the number of frames contributing to the object counts (see count_scale),

\(D\) is the dark current in electrons per hour (see darkcurr),

\(t_{\rm exp}\) is the effective exposure time in seconds (see exptime),

\(V_{\rm rn}\) is the detector readnoise variance (i.e., read-noise squared; see rn_var),

\(V_{\rm proc}\) is added variance from image processing (e.g., bias subtraction; see proc_var), and

\(\epsilon\) is an added error term that imposes a maximum signal-to-noise on the observed counts.

This function consolidates terms that do not change with the forward modeling of the sky + object counts. That is, this function calculates

\[V_{\rm base} = s^2\ \left[ D t_{\rm exp} / 3600 + V_{\rm rn} + V_{\rm proc} \right]\]

such that the first equation can be re-written as

\[V = s {\rm max}(0,c) + V_{\rm base} + \epsilon^2 {\rm max}(0, c)^2.\]

Warning

If \(s\) (count_scale) is provided, the variance will be 0 wherever \(s \leq 0\).
Note that dark current is typically given in electrons per second per pixel. If on-chip binning was used for the detector readout, each binned pixel will have accummulated the expected dark-current (in e-/s/pixel) multiplied by the number of binned pixels. Beware the units of darkcurr, both in that it is dark-current per hour and that it is the dark-current expected in the binned pixel. For example, see the calling function pypeit.images.rawimage.RawImage.build_ivar().
The input arrays can have any dimensionality (i.e., they can be single 2D images or a 3D array containing multiple 2D images); however, the exposure time must be a scalar applied to all array values.

Parameters:

rn_var (numpy.ndarray) – A 2D array with the readnoise variance (i.e., readnoise squared) from the instrument detector; see rn2_frame(). This should include digitization noise and any difference in the readnoise across the detector due to the use of multiple amplifiers. Readnoise should be in e-, meaning this is in elections squared.
darkcurr (float, numpy.ndarray, optional) – Dark current in electrons per hour (as is the convention for the DetectorContainer object) if the exposure time is provided, otherwise in electrons. Note that this is the dark-current in each read pixel, meaning you likely need to multiply the quoted detector dark-current by the number of pixels in a bin (e.g., 4 for 2x2 binning) for binned data. If None, set to 0. If a single float, assumed to be constant across the full image. If an array, the shape must match rn_var.
exptime (float, optional) – Exposure time in seconds. If None, dark current must be in electrons.
proc_var (float, numpy.ndarray, optional) – Additional variance terms to include that are due to the image processing steps (e.g., bias subtraction). If None, set to 0. If a single float, assumed to be constant across the full image. If an array, the shape must match rn_var.
count_scale (float, numpy.ndarray, optional) – A scale factor that has already been applied to the provided counts. It accounts for the number of frames contributing to the provided counts, and the relative throughput factors that can be measured from flat-field frames. For example, if the image has been flat-field corrected, this is the inverse of the flat-field counts. If None, set to 1. If a single float, assumed to be constant across the full image. If an array, the shape must match rn_var. The variance will be 0 wherever \(s \leq 0\), modulo the provided noise_floor.

Returns:

Base-level variance image computed via the equation above with the same shape as rn_var.

Return type:

numpy.ndarray

pypeit.core.procimg.boxcar_average(arr, boxcar)[source]

Boxcar average an array.

The dimensionality of the boxcar matches the dimensionality of the provided array. Averages are only performed over the valid region of the array; see the return description. In contrast with a boxcar smoothing algorithm, the step between each boxcar average is the size of the box, not one pixel.

Parameters:

arr (numpy.ndarray) – Array to average. Can be any shape and dimensionality.
boxcar (int, tuple) – Integer number of pixels to average. If a single integer, all axes are averaged with the same size box. If a tuple, the integer is defined separately for each array axis; length of tuple must match the number of array dimensions.

Returns:

The averaged array. If boxcar is a single integer, the returned array shape is:

tuple([s//boxcar for s in arr.shape])

A similar operation gives the shape when boxcar has elements defined for each array dimension. If the input array is not an integer number of boxcar pixels along a given dimension, the remainder of the array elements along that dimension are ignored (i.e., pixels within the modulus of the array shape and boxcar of the end of the array dimension are ignored).

Return type:

numpy.ndarray

pypeit.core.procimg.boxcar_fill(img, width, bpm=None, maxiter=None, fill_value=nan)[source]

Use convolution with a boxcar kernel to iteratively fill masked regions of an image.

Warning

Depending the size of the kernel, the masked regions, and the maximum number of iterations, the procedure may not be able to fill all masked pixels. Any left-over masked pixels are returned with the provided fill_value.

Parameters:

img (numpy.ndarray, numpy.ma.MaskedArray) – Image to fill. Must be 2D. If an unmasked array, function will simply return img if mask is not provided.
width (int) – Width of the square box to use for the convolution kernel. Must be odd.
bpm (numpy.ndarray, optional) – Image bad-pixel mask. If img is a masked array, this mask is combined with the img mask attribute.
maxiter (int, optional) – The maximum number of fill iterations. If None, iterations continue until all masked pixels are filled.
fill_value (float, optional) – If the number of allowed iterations are insufficient to fill the array, replace masked pixels with this value.

Returns:

Array with the same shape as img with the masked pixels filled by the average of the unmasked pixels in the boxcar kernel.

Return type:

numpy.ndarray

pypeit.core.procimg.boxcar_replicate(arr, boxcar)[source]

Boxcar replicate an array.

Parameters:

arr (numpy.ndarray) – Array to replicate. Can be any shape and dimensionality.
boxcar (int, tuple) – Integer number of times to replicate each pixel. If a single integer, all axes are replicated the same number of times. If a tuple, the integer is defined separately for each array axis; length of tuple must match the number of array dimensions.

Returns:

The block-replicated array.

Return type:

numpy.ndarray

pypeit.core.procimg.cr_screen(a, mask_value=0.0, spatial_axis=1)[source]

Calculate the significance of pixel deviations from the median along the spatial direction.

No type checking is performed of the input array; however, the function assumes floating point values.

Parameters:

a (numpy.ndarray) – Input 2D array
mask_value (float) – (Optional) Values to ignore during the calculation of the median. Default is 0.0.
spatial_axis (int) – (Optional) Axis along which to calculate the median. Default is 1.

Returns:

Returns a map of \(|\Delta_{i,j}|/\sigma_j\), where \(\Delta_{i,j}\) is the difference between the pixel value and the median along axis \(i\) and \(\sigma_j\) is robustly determined using the median absolute deviation, \(sigma_j = 1.4826 MAD\).

Return type:

numpy.ndarray

pypeit.core.procimg.gain_frame(amp_img, gain)[source]

Generate an image with the gain for each pixel.

Parameters:

amp_img (numpy.ndarray) – Integer array that identifies which (1-indexed) amplifier was used to read each pixel.
gain (array-like) – List of amplifier gain values in e-/ADU. Must be that the gain for amplifier 1 is provided by gain[0], etc.

Returns:

Image with the gain for each pixel.

Return type:

numpy.ndarray

pypeit.core.procimg.grow_mask(mask, radius)[source]

Grow pixels flagged as True in a boolean mask by the provided radius.

This is largely a convience wrapper for scipy.ndimage.binary_dilation.

Parameters:

mask (numpy.ndarray) – Boolean mask to process. Pixels flagged as True are expanded into circles with the provided radius.
radius (scalar-like) – Radius in pixels to grow the mask.

Returns:

The boolean mask grown with the masked region grown by the provided radius.

Return type:

numpy.ndarray

pypeit.core.procimg.lacosmic(sciframe, saturation=None, nonlinear=1.0, bpm=None, varframe=None, maxiter=1, grow=1.5, remove_compact_obj=True, sigclip=5.0, sigfrac=0.3, objlim=5.0, rm_false_pos=True)[source]

Identify cosmic rays using the L.A.Cosmic algorithm.

See Peter van Dokkum’s L.A.Cosmic website and van Dokkum (2001, PASP, 113, 1420).

This routine is mostly courtesy of Malte Tewes with some updates/alterations by the PypeIt developers.

Parameters:

sciframe (numpy.ndarray) – Science frame to process.
saturation (float, numpy.ndarray, optional) – The saturation level of the detector in units that match the provided frame. This is used to flag pixels to ignore. Can be a single value or an array; if the latter, the shape must match sciframe. If None, no pixels are flagged as being saturated.
nonlinear (float, numpy.ndarray, optional) – The fraction of the saturation level at which the detector response becomes non-linear. Can be a single value or an array; if the latter, the shape must match sciframe.
bpm (numpy.ndarray, optional) – Bad-pixel mask for the input science frame. If None, all pixels are considered available to be flagged as cosmic rays. Shape must match sciframe.
varframe (numpy.ndarray, optional) – The variance in the science frame to process. If None, the variance is estimated by the absolute value of a \(5\times 5\) boxcar median filter of the image.
maxiter (int, optional) – Maximum number of detection iterations.
grow (float, optional) – The radius (in pixels) used to grow the region masked around detected cosmic rays. See grow_mask(). If \(\leq 0\), the cosmic-ray mask regions are not grown.
remove_compact_obj (bool, optional) – Remove likely compact objects from the set of detected cosmic rays. This is performed by default here and in the original L.A.Cosmic algorithm.
sigclip (float, optional) – Threshold for identifying a cosmic ray
sigfrac (float, optional) – Fraction of sigclip used to define the lower threshold used for pixels neighboring identified cosmic-rays.
objlim (float, optional) – Contrast limit between a cosmic-ray and an underlying object
rm_false_pos (bool, optional) – Apply algorithm to detect and remove false positives. This is not a traditional component of the L.A.Cosmic algorithm.

Returns:

Boolean array flagging pixels with detected cosmic rays; True means the pixel has a cosmic ray.

Return type:

numpy.ndarray

pypeit.core.procimg.pattern_frequency(frame, axis=1)[source]

Using the supplied 2D array, calculate the pattern frequency along the specified axis.

Parameters:

frame (numpy.ndarray) – 2D array to measure the pattern frequency
axis (int, optional) – Which axis should the pattern frequency be measured?

Returns:

The frequency of the sinusoidal pattern.

Return type:

float

pypeit.core.procimg.rect_slice_with_mask(image, mask, mask_val=1)[source]

Generate rectangular slices from a mask image.

Parameters:

image (numpy.ndarray) – Image to mask
mask (numpy.ndarray) – Mask image
mask_val (int, optional) – Value to mask on

Returns:

The image at mask values and a 2-tuple with the slice objects that select the masked data.

Return type:

tuple

pypeit.core.procimg.replace_column_linear(img, left, right)[source]

Replace the column values between left and right indices for all rows by a linear interpolation between the columns just outside the region.

If possible, extrapolation is used for columns at the end of the image with no left or right reference column (left==0 or right==img.shape[1]) using the two most adjacent columns. Otherwise, this function calls replace_column_mean().

Parameters:

img (numpy.ndarray) – Image with values to both use and replace.
left (int) – Inclusive starting column index.
right (int) – Exclusive ending column index.

pypeit.core.procimg.replace_column_mean(img, left, right)[source]

Replace the column values between left and right indices for all rows by the mean of the columns just outside the region.

Columns at the end of the image with no left or right reference column (left==0 or right==img.shape[1]) are just replaced by the closest valid column.

Parameters:

img (numpy.ndarray) – Image with values to both use and replace.
left (int) – Inclusive starting column index.
right (int) – Exclusive ending column index.

pypeit.core.procimg.replace_columns(img, bad_cols, replace_with='mean', copy=False)[source]

Replace bad image columns.

Parameters:

img (numpy.ndarray) – A 2D array with image values to replace.
bad_cols (numpy.ndarray) – Boolean array selecting bad columns in img. Must have the correct shape.
replace_with (str, optional) – Method to use for the replacements. Can be ‘mean’ (see replace_column_mean()) or ‘linear’ (see replace_column_linear()).
copy (bool, optional) – Copy img to a new array before making any modifications. Otherwise, img is modified in-place.

Returns:

The modified image, which is either a new array or points to the in-place modification of img according to the value of copy.

Return type:

numpy.ndarray

pypeit.core.procimg.rn2_frame(datasec_img, ronoise, units='e-', gain=None, digitization=False)[source]

Construct a readnoise variance image.

Provided the detector readnoise and gain for each amplifier, this constructs an image with the combination of the readnoise and digitization (or quantization) noise expected for a single detector readout. Digitization noise is a fixed \(\sqrt{1/12}\) ADU [1] [2], derived as the second moment of a uniform distribution between values of -1/2 to 1/2 (i.e., the variance associated with converting a number of electrons into an ADU integer quantized by the gain). The digitization noise is typically much smaller than the readnoise, unless the gain is very large, and, depending on how it was measured, the digitization noise is most often incorporated in the documented readnoise of the given instrument. To include the digitization noise in the variance, you must provide gain and set digitization=True.

The variance calculation in electrons is \(V = {\rm RN}^2 + \gamma^2/12\), when including the digitization noise, and simply \(V = {\rm RN}^2\) otherwise; where RN is the readnoise and \(\gamma\) is the gain in e-/ADU. In the rare case one would need the units in ADU, the returned variance is \(V/\gamma^2\).

Parameters:

datasec_img (numpy.ndarray) – An integer array indicating the 1-indexed amplifier used to read each pixel in the main data section of the detector. Values of 0 are ignored. Amplifier numbers are expected sequential and match the number of readnoise and gain values provided. The shape of this image dictates the shape of the output readnoise variance image.
ronoise (float, array-like) – The value of the readnoise for each amplifier in electrons (e-). If there is only one amplifier, this can be provided as a single float.
units (str, optional) – Units for the output variance. Options are 'e-' for variance in square electrons (counts) or 'ADU' for square ADU.
gain (float, array-like, optional) – The value of the gain for each amplifier in e-/ADU. If digitization is False, this is ignored.
digitization (bool, optional) – Include digitization error in the calculation. If True, gain must be provided.

Returns:

The image variance resulting from reading the detector in the selected units for each pixel. The shape is the same as datasec_img. Pixels where datasec_img is 0 are set to 0.

Return type:

numpy.ndarray

pypeit.core.procimg.subtract_overscan(rawframe, datasec_img, oscansec_img, method='savgol', params=[5, 65], var=None)[source]

Subtract the overscan

Possible values of method:

polynomial: Fit a polynomial to the overscan region and subtract it.
savgol: Use a Savitzky-Golay filter to fit the overscan region and
subtract it.
median: Use the median of the overscan region to subtract it.
odd_even: Use the median of the odd and even rows/columns to subtract (MDM/OSMOS)

Parameters:

rawframe (numpy.ndarray) – Frame from which to subtract overscan. Must be 2d.
datasec_img (numpy.ndarray) – An array the same shape as rawframe that identifies the pixels associated with the data on each amplifier; 0 for no data, 1 for amplifier 1, 2 for amplifier 2, etc.
oscansec_img (numpy.ndarray) – An array the same shape as rawframe that identifies the pixels associated with the overscan region on each amplifier; 0 for no data, 1 for amplifier 1, 2 for amplifier 2, etc.
method (str, optional) – The method used to fit the overscan region. Options are chebyshev, polynomial, savgol, median. (“polynomial” is deprecated and will be removed)
params (list, optional) – Parameters for the overscan subtraction. For method=chebyshev or method=polynomial``set ``params to the order; for method=savgol, set params to the order and window size; for method=median, params are ignored.
var (numpy.ndarray, optional) – Variance in the raw frame. If provided, must have the same shape as rawframe and used to estimate the error in the overscan subtraction. The estimated error is the standard error in the median for the pixels included in the overscan correction. This estimate is also used for the 'savgol' method as an upper limit. If None, no variance in the overscan subtraction is calculated, and the 2nd object in the returned tuple is None.

Returns:

The input frame with the overscan region subtracted and an estimate of the variance in the overscan subtraction; both have the same shape as the input rawframe. If var is not provided, the 2nd returned object is None.

Return type:

tuple

pypeit.core.procimg.subtract_pattern(rawframe, datasec_img, oscansec_img, frequency=None, axis=1)[source]

Subtract a sinusoidal pattern from the input rawframe. The algorithm calculates the frequency of the signal, generates a model, and subtracts this signal from the data. This sinusoidal pattern noise was first identified in KCWI, but the source of this pattern noise is not currently known. The pattern model is generated with a three step process:

Given a first guess at the frequency, calculate how frequency depends on pixel row (slight linear dependence)

Using the frequency model, calculate the amplitude of the signal (usually a constant for all pixel rows)

Using the model of the frequency and the amplitude, calculate the phase of the signal for each pixel row. Note that the phase is different for each pixel row.

A complete detector model is generated for the pattern noise using the frequency+amplitude+phase, and an estimate of the improved effective read noise is provided.

Parameters:

rawframe (numpy.ndarray) – Frame from which to subtract overscan
numamplifiers (int) – Number of amplifiers for this detector.
datasec_img (numpy.ndarray) – An array the same shape as rawframe that identifies the pixels associated with the data on each amplifier. 0 for not data, 1 for amplifier 1, 2 for amplifier 2, etc.
oscansec_img (numpy.ndarray) – An array the same shape as rawframe that identifies the pixels associated with the overscan region on each amplifier. 0 for not data, 1 for amplifier 1, 2 for amplifier 2, etc.
frequency (float, list, optional) – The frequency (or list of frequencies - one for each amplifier) of the sinusoidal pattern. If None, the frequency of each amplifier will be determined from the overscan region.
axis (int, optional) – Which axis should the pattern subtraction be applied?

Returns:

The input frame with the pattern subtracted

Return type:

numpy.ndarray

pypeit.core.procimg.trim_frame(frame, mask)[source]

Trim the masked regions from a frame.

Parameters:

frame (numpy.ndarray) – Image to be trimmed
mask (numpy.ndarray) – Boolean image set to True for values that should be trimmed and False for values to be returned in the output trimmed image.

Returns:

Trimmed image

Return type:

numpy.ndarray

Raises:

PypeItError – Error raised if the trimmed image includes masked values because the shape of the valid region is odd.

pypeit.core.procimg.variance_model(base, counts=None, count_scale=None, noise_floor=None)[source]

Calculate the expected variance in an image.

The full variance model (see variance_model()), \(V\), is:

\[V = s^2\ \left[ {\rm max}(0, C) + D t_{\rm exp} / 3600 + V_{\rm rn} + V_{\rm proc} \right] + \epsilon^2 {\rm max}(0, c)^2\]

where:

\(c=s\ C\) are the rescaled observed sky + object counts (see counts),

\(C\) is the observed number of sky + object counts,

\(s=s\prime / N_{\rm frames}\) is a scale factor derived from the (inverse of the) flat-field frames plus the number of frames contributing to the object counts (see count_scale),

\(D\) is the dark current in electrons per hour,

\(t_{\rm exp}\) is the effective exposure time in seconds,

\(V_{\rm rn}\) is the detector readnoise variance (i.e., read-noise squared),

\(V_{\rm proc}\) is added variance from image processing (e.g., bias subtraction), and

\(\epsilon\) is an added error term that imposes a maximum signal-to-noise on the observed counts (see noise_floor).

The function base_variance() consolidates all terms that do not change with the forward modeling of the sky + object counts into a single “base-level” variance

\[V_{\rm base} = s^2\ \left[ D t_{\rm exp} / 3600 + V_{\rm rn} + V_{\rm proc} \right]\]

such that the first equation can be re-written as

\[V = s {\rm max}(0,c) + V_{\rm base} + \epsilon^2 {\rm max}(0, c)^2,\]

which is the quantity returned by this function.

We emphasize that this is a model for the per-pixel image variance. In real data, the as-observed pixel values are used to estimate the Poisson error in the observed counts. Because of the variance in the image, this systematically overestimates the variance toward low counts (\(\lesssim 2 \sigma_{\rm rn}\)), with a bias of approximately \(1.4/\sigma_{\rm rn}\) for \(C=0\) (i.e., about 20% for a readnoise of 2 e-) and less than 10% for \(C=1\).

Warning

If \(s\) (count_scale) is provided, the variance will be 0 wherever \(s \leq 0\), modulo the provided noise_floor.
The input arrays can have any dimensionality (i.e., they can be single 2D images or a 3D array containing multiple 2D images); however, currently the noise floor must be a single scalar applied to all array values.

Parameters:

base (numpy.ndarray) – The “base-level” variance in the data set by the detector properties and the image processing steps. See base_variance(); \(V_{\rm base}\) in the equations above.
counts (numpy.ndarray, optional) – An array with the number of source-plus-sky counts, possibly rescaled by a relative throughput; see \(c\) in the equations above. Because this is used to calculate the noise floor, this must be provided if noise_floor is not None. Shape must match base. Shape must match base.
count_scale (float, numpy.ndarray, optional) – A scale factor that has already been applied to the provided counts; see \(s\) in the equations above. It accounts for the number of frames contributing to the provided counts, and the relative throughput factors that can be measured from flat-field frames. For example, if the image has been flat-field corrected, this is the inverse of the flat-field counts. If None, no scaling is expected, meaning counts are exactly the observed detector counts. If a single float, assumed to be constant across the full image. If an array, the shape must match base. The variance will be 0 wherever \(s \leq 0\), modulo the provided noise_floor.
noise_floor (float, optional) – A fraction of the counts to add to the variance, which has the effect of ensuring that the S/N is never greater than 1/noise_floor; see \(epsilon\) in the equations above. If None, no noise floor is added. If not None, counts must be provided.

Returns:

Variance image computed via the equation above with the same shape as base.

Return type:

numpy.ndarray