JP-3622 Update refpix step for NIR detectors to use convolution kernel (

#8726)

Co-authored-by: David Law <[email protected]>
Co-authored-by: Melanie Clarke <[email protected]>
Co-authored-by: Tyler Pauly <[email protected]>

changes/8726.refpix.rst

@@ -0,0 +1 @@
+Implemented SIRS algorithm instead of running median for side pixels of NIR full-frame data. Running median is still default.

docs/jwst/refpix/arguments.rst

@@ -52,3 +52,25 @@ in IRS2 mode will be processed along with the normal pixels and preserved
 in the output.  This option is intended for calibration or diagnostic reductions
 only. For normal science operation, this argument should always be False,
 so that interleaved pixels are stripped before continuing processing.
+
+*  ``--refpix_algorithm``
+
+The ``refpix_algorithm`` argument is only relevant for all NIR full-frame
+data, and can be set to 'median' (default) to use the running median or
+'sirs' to use the Simple Improved Reference Subtraction (SIRS).
+
+*  ``--sigreject``
+
+The ``sigreject`` argument is the number of sigmas to reject as outliers in the
+SIRS algorithm. The value is expected to be a float.
+
+*  ``--gaussmooth``
+
+The ``gaussmooth`` argument is the width of Gaussian smoothing kernel to use as
+a low-pass filter. The numerical value is expected to be a float.
+
+*  ``--halfwidth``
+
+The ``halfwidth`` argument is the half-width of convolution kernel to build. The
+numerical value is expected to be an integer.
+

docs/jwst/refpix/description.rst

@@ -76,6 +76,18 @@ NIR Detector Data
        (set by the step parameter ``side_gain``, which defaults to 1.0) is
        subtracted from the full group on a row-by-row basis.  Note that the ``odd_even_rows``
        parameter is ignored for NIR data when the side reference pixels are processed.
+       If the ``--refpix_algorithm`` option is set to 'sirs', the Simple Improved
+       Reference Subtraction (SIRS) method will be used instead of the running median.
+       The SIRS revision uses the left and right side reference pixels as described
+       in https://doi.org/10.1117/1.JATIS.8.2.028002. This implementation uses a
+       mathematically equivalent formulation using convolution kernels rather than
+       Fourier transforms, with the convolution kernel truncated where the weights
+       approach zero. There are two convolution kernels for each readout channel,
+       one for each of the left and right reference pixels. These kernels are the
+       Fourier transforms of the weight coefficients (further description in the paper).
+       The approach implemented here makes nearly optimal use of the reference pixels.
+       It reduces read noise by 5-10% relative to a running median filter of the
+       reference pixels, and reduces 1/f "striping" noise by a larger factor than this.
     #. Transform the data back to the JWST focal plane, or DMS, frame.
 
 MIR Detector Data

jwst/refpix/optimized_convolution.py

@@ -0,0 +1,180 @@
+#
+# Module for using the Simple Improved Reference Subtraction (SIRS) algorithm
+# to improve the 1/f noise, only for full frame non-IRS2 NIR data
+#
+
+import logging
+import numpy as np
+
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
+
+
+def make_kernels(sirs_kernel_model, detector, gaussmooth, halfwidth):
+    """
+    Make convolution kernels from Fourier coefficients in the reference file.
+
+    Parameters:
+    -----------
+
+    sirs_kernel_model : `~jwst.datamodels.SIRSKernelModel`
+        Data model containing the Fourier coefficients from the reference files for
+        Simple Improved Reference Subtraction (SIRS)
+
+    detector : str
+        Name of the detector of the input data
+
+    gaussmooth : float
+        Width of Gaussian smoothing kernel to use as a low-pass filter on reference file's coefficients
+
+    halfwidth : int
+        Half-width of convolution kernel to build from reference file's coefficients
+
+    Returns:
+    --------
+    kernels: list
+        List of kernels appropriate for convolution with the left and right reference pixels.
+
+    """
+
+    gamma, zeta = get_conv_kernel_coeffs(sirs_kernel_model, detector)
+    if gamma is None or zeta is None:
+        log.info(f'Optimized convolution kernel coefficients NOT found for detector {detector}')
+        return None
+
+    kernels_left = []
+    kernels_right = []
+    for chan in range(gamma.shape[0]):
+        n = len(gamma[chan]) - 1
+        kernel_left = np.fft.fftshift(np.fft.irfft(gamma[chan]))[n - halfwidth:n + halfwidth + 1]
+        kernel_right = np.fft.fftshift(np.fft.irfft(zeta[chan]))[n - halfwidth:n + halfwidth + 1]
+
+        x = np.arange(-halfwidth, halfwidth + 1)
+        window = np.exp(-x ** 2 / (2 * gaussmooth ** 2))
+        window /= np.sum(window)
+
+        kernel_right = np.convolve(kernel_right, window, mode='same')
+        kernel_left = np.convolve(kernel_left, window, mode='same')
+
+        kernels_right += [kernel_right]
+        kernels_left += [kernel_left]
+
+    return [kernels_left, kernels_right]
+
+
+def get_conv_kernel_coeffs(sirs_kernel_model, detector):
+    """
+    Get the convolution kernels coefficients from the reference file
+
+    Parameters:
+    -----------
+
+    sirs_kernel_model : `~jwst.datamodels.SIRSKernelModel`
+        Data model containing the Fourier coefficients from the reference files for
+        Simple Improved Reference Subtraction (SIRS)
+
+    detector : str
+        Name of the detector of the input data
+
+    Returns:
+    --------
+
+    gamma: numpy array
+        Fourier coefficients
+
+    zeta: numpy array
+        Fourier coefficients
+    """
+    mdl_dict = sirs_kernel_model.to_flat_dict()
+    gamma, zeta = None, None
+    for item in mdl_dict:
+        det = item.split(sep='.')[0]
+        if detector.lower() == det.lower():
+            arr_name = item.split(sep='.')[1]
+            if arr_name == 'gamma':
+                gamma = np.array(mdl_dict[item])
+            elif arr_name == 'zeta':
+                zeta = np.array(mdl_dict[item])
+        if gamma is not None and zeta is not None:
+            break
+    return gamma, zeta
+
+
+def apply_conv_kernel(data, kernels, sigreject=4.0):
+    """
+    Apply the convolution kernel.
+
+    Parameters:
+    -----------
+
+    data : 2-D numpy array
+        Data to be corrected
+
+    kernels : list
+        List containing the left and right kernels
+
+    sigreject: float
+        Number of sigmas to reject as outliers
+
+    Returns:
+    --------
+
+    data : 2-D numpy array
+        Data model with convolution
+    """
+    data = data.astype(float)
+    npix = data.shape[-1]
+
+    kernels_l, kernels_r = kernels
+    nchan = len(kernels_l)
+
+    L = data[:, :4]
+    R = data[:, -4:]
+
+    # Find the approximate standard deviations of the reference pixels
+    # using an outlier-robust median approach. Mask pixels that differ
+    # by more than sigreject sigma from this level.
+    # NOTE: The Median Absolute Deviation (MAD) is calculated as the
+    #   median of the absolute differences between data values and their
+    #   median. For normal distribution MAD is equal to 1.48 times the
+    #   standard deviation but is a more robust estimate of the dispersion
+    #   of data values.The calculation of MAD is straightforward but
+    #   time-consuming, especially if MAD estimates are needed for the
+    #   local environment around every pixel of a large image. The
+    #   calculation is MAD = np.median(np.abs(x-np.median(x))).
+    #   Reference: https://www.interstellarmedium.org/numerical_tools/mad/
+    MAD = 1.48
+    medL = np.median(L)
+    sigL = MAD * np.median(np.abs(L - medL))
+    medR = np.median(R)
+    sigR = MAD * np.median(np.abs(R - medR))
+
+    # nL and nR are the number of good reference pixels in the left and right
+    # channel in each row. These will be used in lieu of replacing the values
+    # of those pixels directly.
+    goodL = 1 * (np.abs(L - medL) <= sigreject * sigL)
+    nL = np.sum(goodL, axis=1)
+    goodR = 1 * (np.abs(R - medR) <= sigreject * sigR)
+    nR = np.sum(goodR, axis=1)
+
+    # Average of the left and right channels, replacing masked pixels with zeros.
+    # Appropriate normalization factors will be computed later.
+    L = np.sum(L * goodL, axis=1) / 4
+    R = np.sum(R * goodR, axis=1) / 4
+    for chan in range(nchan):
+        kernel_l = kernels_l[chan]
+        kernel_r = kernels_r[chan]
+
+        # Compute normalizations so that we don't have to directly
+        # replace the values of flagged/masked reference pixels.
+        normL = np.convolve(np.ones(nL.shape), kernel_l, mode='same')
+        normL /= np.convolve(nL / 4, kernel_l, mode='same')
+        normR = np.convolve(np.ones(nR.shape), kernel_r, mode='same')
+        normR /= np.convolve(nR / 4, kernel_r, mode='same')
+        template = np.convolve(L, kernel_l, mode='same') * normL
+        template += np.convolve(R, kernel_r, mode='same') * normR
+        data[:, chan * npix * 3 // 4:(chan + 1) * npix * 3 // 4] -= template[:, np.newaxis]
+
+    log.debug('Optimized convolution kernel applied')
+    return data
+