pypeit.core.moment module

Module to compute moments.

pypeit.core.moment.moment1d(flux, col, width, ivar=None, bpm=None, fwgt=None, row=None, weighting='uniform', order=0, bounds=None, fill_error=-1.0, mesh=False)

Compute one-dimensional moments of the provided image within an aperture along its second axis (axis=1).

This method allows for computations of the zeroth, first, and second moments (see order). The aperture used for the calculation is centered at the provided col pixel with a width defined by width; however, this definition depends on the type of weighting applied (see weighting). Formulae for each moment are as follows. The zeroth moment (order=0) computes the discrete weighted sum of the flux within the aperture:

\[\begin{split}\mu_0 &= \sum_i w_i f_i \\ \epsilon_{\mu_0}^2 &= \sum_i (w_i \epsilon_{f,i})^2,\end{split}\]

where \(f\) is the flux in each pixel \(i\), \(\epsilon_f\) is its error, and \(w\) is the assigned weight (see weighting). The first moment (order=1) computes the flux-weighted center of window:

\[\begin{split}\mu_1 &= \frac{\sum_i x_i w_i f_i }{\mu_0} \\ \epsilon_{\mu_1}^2 &= \mu_0^{-2}\ \sum_i [ w_i \epsilon_{f,i}\ (x_i - \mu_1)]^2,\end{split}\]

where \(x\) is the pixel position along the 2nd axis (see col). The second moment (order=2) computes the variance of the flux profile about its center within the window:

\[\begin{split}\mu_2 &= \frac{\sum_i x^2_i w_i f_i }{\mu_0} - \mu_1^2 \\ \epsilon_{\mu_2}^2 &= \mu_0^{-2}\ \sum_i w^2_i \epsilon^2_{f,i}\ [(x_i - \mu_1)^2 + \mu_2]^2.\end{split}\]

The values returned for the second-moment calculation are actually the standard deviation instead of the variance, where:

\[\begin{split}\sigma = \sqrt{\mu_2} \\ \epsilon_{\sigma} = \frac{\epsilon_{\mu_2}}{2 \sigma}.\end{split}\]

The method uses numpy.ma.MaskedArray objects to keep track of math errors, such as divisions by 0. The returned boolean array indicates when these errors occurred, and the method replaces these errors with the original centers.

The shape of arrays depend on the shapes of the input image (flux), columns about which to perform the calculation (col), and rows from which to select the column data (row) as listed below. In the following, the input flux array must be 2D and has shape \((N_{\rm row}, N_{\rm col})\), whereas the number of rows and columns identified for a moment calculation are \(N_{\rm mom,row}\) and \(N_{\rm mom,col}\), respectively. The possible shapes of the output arrays are:

• If row is None and col is 1D, the moments are calculated centered at each col for all rows in flux. The shape of the array for each moment is then \((N_{\rm row}, N_{\rm mom,col})\)

• If row is None and col is 2D, col must have shape \((N_{\rm row}, N_{\rm mom,col})\) and the shape of the output array per moment calculated is the same.

• If row is a single integer and col is 1D, the moments are calculated for each col only at the provided row in flux. The shape of the array for each moment is then \((N_{\rm mom,col},)\).

• If row is a 1D integer array and col is 1D, the default behavior is for the moments to be calculated at each col for each provided row. The shape of the array for each moment would then be \((N_{\rm mom,row},N_{\rm mom,col})\). This will be true even if \(N_{\rm mom,row} = N_{\rm mom,col}\), as long as mesh is True. However, if \(N_{\rm mom,row} = N_{\rm mom,col}\) and mesh is False, only the matched pairs of col and row will be used and the output shape per moment will be \((N_{\rm mom,row},)\).

• If row is a 1D integer array and col is 2D, the first axis of col must be the same length as row. The shape of the array for each moment is then \((N_{\rm mom,row},N_{\rm mom,col})\), which is the same as the input col.

• If row is 2D and col is 2D, they must have the same shape.

• If row is a single integer and col is 2D, or if row is 2D and col is 1D, the method raises an error.

Note

• This is an entirely general function, as reflected by the nomenclature used in the call. As used within PypeIt, the PypeIt image orientation convention means that moments are always taken along the spatial direction; i.e., col is the spatial coordinate and row is the spectral coordinate.

• This function is a generalization of and builds on the heritage of functions in idlspec2d, specifically trace_fweight, trace_gweight, extrace_asymbox2, extract_boxcar.

Warning

The function has significant setup/input checking. Most of this introduces limited overhead with the exception of the handling of ivar, bpm, and fwgt. If any of these are provided as None on input, an array is constructed (unity for ivar and fwgt and all False for bpm) that serves as a place holder. If repetitive calls to the function are expected and any of these arrays are missing, significant efficiency gains can be made by providing pre-built values for these arrays so that time isn’t lost in allocating the placeholder arrays in every call.

Todo

Optimize the code for efficiency, regardless of the input.

Parameters:

• flux (numpy.ndarray) – Intensity image with shape \((N_{\rm row}, N_{\rm col})\).

• col (numpy.ndarray) – Floating-point center along the 2nd axis for the integration window in pixel index (first value located at index 0). This can either be a 1D or 2D array. See restrictions on the shape in the description above.

• width (float, numpy.ndarray) – The meaning of the parameter depends on the value of weighting. If weighting==’uniform’, the width of the integration window in columns, centered at the input col. If weighting==’gaussian’, the \(\sigma\) of a pixelated Gaussian weighting function. The width of the integration window in columns, centered at the input col, is always 6*width (i.e., the half-width is \(3\sigma\)). The provided value can be a scalar to use a constant window definition, or it can be an array variable integration window where the array must have the same shape as col.

• ivar (numpy.ndarray, optional) – Inverse variance of the image intensity. If not provided, unity variance is used. If provided, must have the same shape as flux.

• bpm (numpy.ndarray, optional) – Boolean bad-pixel mask for the input image. True values are ignored, False values are included. If not provided, all pixels are included. If provided, must have the same shape as flux.

• fwgt (numpy.ndarray, optional) – An additional weight to apply to each pixel in flux. If None, weights are uniform. Otherwise, the \(w_i\) from above are the product of this weight and the result of the scheme set using the weighting argument.

• row (int, numpy.ndarray, optional) – Integer or integer array with the position along the first axis (axis=0) for the moment calculation. This can either be None, an integer, or a 1D or 2D array. See restrictions on the shape in the description above.

• weighting (str, optional) – The weighting to apply to the position within each integration window (see width above). This must be (case-insensitive) either ‘uniform’ for uniform weighting or ‘gaussian’ for weighting by a Gaussian centered at the input guess coordinates and integrated over the pixel width.

• order (int, array-like, optional) – The order of the moment(s) to calculate. Can be a single integer or a list. Moments to calculate must be 0, 1, or 2; at most order can be [0,1,2]. The shape of the output arrays depends on the number of moments calculated. Note that the calculation of the orders is necessarily sequential; i.e., setting order=2 means that the zeroth and first moments have to be calculated anyway. The order must be provided in sorted order; i.e., you cannot pass order=[2,1].

• bounds (tuple, optional) –

  A two-tuple with the lower and upper limit for each moment order. If None, no bounds are imposed. If not None, an upper and lower bound must be provided for each moment to compute; i.e., if more than one moment is computed, each element of the two-tuple must be an array-like object that matches the length of order. To set an upper or lower bound only, set the unbounded component to None. Bounds for the zeroth and second order moments are in an absolute sense, whereas first-order bounds are relative to the input col. Measurements that hit the bounds are masked; see the description of the returned objects. For example, to flag anything without a positive zeroth moment or a maximum shift from the input center of 1 pixel, call the method with arguments:

  order=[0,1], bounds=([0,-1], [None,1])
  

• fill_error (float, optional) – Value to use as filler for undetermined moments, resulting from either the input bad-pixel mask or computational issues (division by zero, etc.; see return description below).

• mesh (bool, optional) – If col and row are 1D vectors of the same length, this determines if each col and row should be paired (mesh is False) or used to construct a grid, i.e., every col combined with every row (mesh is True). See the method description.

Returns:

Three numpy.ndarray objects are returned. If more than one moment order is requested, the moments are ordered along the first axis; e.g., if order=[0,1] the outputs moment[0] and moment[1] contain the zeroth and first moments, respectively. The subsequent dimensions of the output arrays are dictated by the input row and col; see the method description. The returned arrays are:

• The moment calculated along the 2nd axis of the input image (axis=1). Masked values (indicated by the third object returned) are 0 for the zeroth and second moments and equal to the input col value for the first moments.

• The formal propagated error (see equations above) in the moments. Errors are only meaningful if ivar is provided. Masked values (indicated by the third object returned) are set to fill_error.

• A boolean bad-pixel mask for output data; True values should be ignored, False values are valid measurements.

Raises:

ValueError – Raised if input shapes are not correct or if the selected weighting is unknown. See method description.

Examples

First setup an image with some Gaussians:

>>> from pypeit.core.moment import moment1d
>>> import numpy as np
>>> from scipy.special import erf
>>> def gauss_comb():
...     c = [45,50,55]
...     img = np.zeros((len(c),100), dtype=float)
...     x = np.arange(100)
...     sig = 5.
...     for i,_c in enumerate(c):
...         img[i,:] = (erf((x-c[i]+0.5)/np.sqrt(2)/sig)
...                         - erf((x-c[i]-0.5)/np.sqrt(2)/sig))/2.
...     return img
...
>>> img = gauss_comb()

Calculate all moments at one column and row:

>>> mu, mue, flag = moment1d(img, 50, 40., row=0, order=[0,1,2])
>>> print(mu)
[ 0.99858297 45.02314924  4.97367636]

Calculate all moments at one column for all rows:

>>> moment1d(img, 50, 40., order=[0,1,2])[0]
array([[ 0.99858297,  0.99993125,  0.99858297],
       [45.02314924, 50.        , 54.97685076],
       [ 4.97367636,  5.00545947,  4.97367636]])

Calculate zeroth moments in all rows centered at column 50

>>> moment1d(img, 50, 40., order=0)[0]
array([0.99858297, 0.99993125, 0.99858297])

Calculate zeroth moments in all rows for three column positions

>>> moment1d(img, [45,50,55], 40., order=0, mesh=True)[0]
array([[0.99993125, 0.99858297, 0.97670951],
       [0.99858297, 0.99993125, 0.99858297],
       [0.97670951, 0.99858297, 0.99993125]])

Calculate first moments in all rows for the three column positions

>>> moment1d(img, [45,50,55], 40., order=1, mesh=True)[0]
array([[45.        , 45.02314924, 45.2814655 ],
       [49.97685076, 50.        , 50.02314924],
       [54.7185345 , 54.97685076, 55.        ]])

Or pick the same column for all rows

>>> moment1d(img, 50, 40., row=[0,1,2], order=1)[0]
array([45.02314924, 50.        , 54.97685076])

Or pick a column unique to each row

>>> moment1d(img, [43,52,57], 40., row=[0,1,2], order=1)[0]
array([44.99688181, 50.00311819, 55.00311819])