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 functionpypeit.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 theDetectorContainer
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 matchrn_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 matchrn_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 matchrn_var
. The variance will be 0 wherever \(s \leq 0\), modulo the providednoise_floor
.
- Returns:
Base-level variance image computed via the equation above with the same shape as
rn_var
.- Return type:
- 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 atuple
, 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:
- 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
ifmask
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 theimg
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:
- 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 atuple
, 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:
- 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:
- 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:
- 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:
- 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 matchsciframe
. 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 matchsciframe
.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. Seegrow_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 raysigfrac (
float
, optional) – Fraction ofsigclip
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 objectrm_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:
- 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:
- 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:
- 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’ (seereplace_column_mean()
) or ‘linear’ (seereplace_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:
- 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 setdigitization=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. Ifdigitization
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 wheredatasec_img
is 0 are set to 0.- Return type:
- 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 asrawframe
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. Formethod=chebyshev
ormethod=polynomial``set ``params
to the order; formethod=savgol
, setparams
to the order and window size; formethod=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
. Ifvar
is not provided, the 2nd returned object is None.- Return type:
- Possible values of
- 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:
- 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:
- 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 providednoise_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 matchbase
. Shape must matchbase
.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, meaningcounts
are exactly the observed detector counts. If a single float, assumed to be constant across the full image. If an array, the shape must matchbase
. The variance will be 0 wherever \(s \leq 0\), modulo the providednoise_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 than1/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: