NISP-S Scaling

I recently realized that the pixel value scaling, implemented in my inSpector code, could be useful more generally, so I have made it available, along with a few other helpful functions, here as a stand-alone Python module. The scaling function scale() was specifically crafted to make most NISP-S (spectral) images look decent without any manual manipulations, as would be required with a FITS image viewer like SAOImage DS9. The code adjusts the histogram of pixel values in such a way that faint regions and very bright regions are simultaneously displayed well. The scaling technique is not exactly histogram equalization nor gamma correction; it uses a combination of arctangent compression, histogram re-shaping, and clipping to achieve the final result. The scale() function can also accept a second, auxiliary image to be scaled based on the parameters computed from another image so that two detector images can be compared consistently.

The histograms below show the original distribution of pixel values in a detector and the output of the scale() function.

Before scaling:

After scaling:

An example of an image that has been scaled using the function is shown below. Before scaling, these images are typically black with just a few visible bright pixels.

Other helpful functions in the module include display() which can be used to display a numpy array containing data from a detector and make_jpegs(), which can create JPEG images of all 16 detectors in an exposure.

Example Usage

The following assumes that the user is in an IPython session and the nispview.py file is either installed or located in the current directory.

import nispview

f = nispview.fits.open('EUC_SIR_W-SCIFRM_0-3_20220930T182125.058129Z.fits')

# open an interactive matplotlib plot showing the scaled image:

nispview.display(f[1].data)

# create an image of each detector in the file:

nispview.make_jpegs('EUC_SIR_W-SCIFRM_0-3_20220930T182125.058129Z.fits')

# Output:

Saving 25463-dither-3-DET11.SCI.jpg.
Saving 25463-dither-3-DET21.SCI.jpg.
Saving 25463-dither-3-DET31.SCI.jpg.
Saving 25463-dither-3-DET41.SCI.jpg.
Saving 25463-dither-3-DET12.SCI.jpg.
Saving 25463-dither-3-DET22.SCI.jpg.
Saving 25463-dither-3-DET32.SCI.jpg.
Saving 25463-dither-3-DET42.SCI.jpg.
Saving 25463-dither-3-DET13.SCI.jpg.
Saving 25463-dither-3-DET23.SCI.jpg.
Saving 25463-dither-3-DET33.SCI.jpg.
Saving 25463-dither-3-DET43.SCI.jpg.
Saving 25463-dither-3-DET14.SCI.jpg.
Saving 25463-dither-3-DET24.SCI.jpg.
Saving 25463-dither-3-DET34.SCI.jpg.
Saving 25463-dither-3-DET44.SCI.jpg.

The acual scale() function looks like this:

def scale(im: np.ndarray, gamma: float=3, hist_peak: float=40, aux_im=None):
    """
    Scales the pixel values in the image array, `im`, such that details in 
    faint regions and bright regions are both clearly visible. Outputs an 
    image with values from 0 and 255, inclusive. Optionally, an auxiliary image can be provided,
    in which case this image will be scaled using the same parameters as the
    primary image, im.

    Returns: Tuple[np.ndarray, Union[np.ndarray, None]]

    One or two scaled images (two if an auxiliary image, `aux_im`, is provided).
    """
    sigmoid = np.arctan
    im = remove_nans(im)
    minval = np.percentile(im, 0.05)
    maxval = np.percentile(im, 99.85)
    data = (im - minval) / (maxval - minval)
    data = 255 * sigmoid(gamma * data) / 1.57079
    left = np.percentile(data, 2)
    vals, bins = np.histogram(data.flat, bins=512, range=(0, 255))
    peak = bins[np.argmax(vals)] - left
    peak = 0.3 if peak <= 0.15 else peak
    correction = interp1d([0,                      peak,             255],
                          [0.08 * hist_peak / peak, hist_peak / peak, 1.0  ],
                          kind='slinear', fill_value=(0,1), bounds_error=False)
    data -= left
    data *= correction(data)
    data = np.clip(data, 0, 255)

    if aux_im is not None:
        aux_data = (aux_im - minval) / (maxval - minval)
        aux_data = 255 * sigmoid(gamma * aux_data) / 1.57079
        aux_data -= left
        aux_data *= correction(aux_data)
        data = np.clip(aux_data, 0, 255)
    else:
        aux_data = None
    return data, aux_data