How to better process and stack JWST data using Python for prettier space pictures?

Question

I thought it would be fun to try my hand at creating beautiful looking astronomy pictures using James Webb Space Telescope (JWST) data which is freely available to download.. if you can make sense of it.

There are tutorials aplenty explaining portions of the task, and crumbs left by others on the same journey which have carried me farther than most I like to think. However I feel I've exhausted what guidance I'm able to find by searching and I'm just not happy with the results.

As a benchmark, I am attempting to recreate the famous Southern Ring Nebula image: https://webbtelescope.org/contents/media/images/2022/033/01G70BGTSYBHS69T7K3N3ASSEB

This seemed a good approach as it's NASA so they kindly explain part of their process.

From that page: "These images are a composite of separate exposures acquired by the James Webb Space Telescope using the NIRCam instrument. Several filters were used to sample narrow and broad wavelength ranges. The color results from assigning different hues (colors) to each monochromatic (grayscale) image associated with an individual filter." This image is a combination of six wavelength exposures.

So I have the exposure data (trying not to call them images at this point) in fits format.

Load, process, display a single exposure

As a first pass at processing and viewing a single exposure:

1. Small border crop since the sensor data doesn't completely fill the array
2. Replace 0s with nans then filter them out which fixes holes in the center of stars
3. Do a very basic brightness mapping since the sensor data doesn't translate directly to a viewable image

from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np
import sys
from astropy.convolution import Gaussian2DKernel, interpolate_replace_nans

fits_filename = sys.argv[1]

hdulist = fits.open(fits_filename, mode="readonly")
h_ref = hdulist[1].header
img_array = hdulist[1].data[80:-80, 80:-80]  # crop

img_array[img_array == 0] = np.nan  # 0s -> nans
kernel = Gaussian2DKernel(x_stddev=3)
img_array = interpolate_replace_nans(img_array, kernel)  # filter out nans

# brightness mapping
xp = [np.nanpercentile(img_array, 4), np.nanpercentile(img_array, 99.9)]
yp = [0, 255]
img_array = np.interp(img_array, xp, yp).astype("uint8")

plt.figure(1, figsize=(6, 6))
im = plt.imshow(img_array, cmap='gray')
plt.show()

Which for the F090W exposure looks like:

"Star alignment" / projection

Next, to stack multiple exposures they need to be aligned since the space telescope is always moving. I've seen "star aligning" techniques, but the sky-coordinate info is in the fits header, so I use reproject to project and interpolate a second fits onto a first reference fits which works well:

from astropy.io import fits
from reproject import reproject_interp

hdulist = fits.open(reference_fits, mode="readonly")
reference_header = hdulist[1].header

hdulist = fits.open(second_fits, mode="readonly")
img_array, footprint = reproject_interp(hdulist[1], reference_header)
img_array = img_array[80:-80, 80:-80]

So of course I built a GUI to load multiple exposures, align, adjust the brightness mappings, assign colors for each "layer" and then stack them to produce a color image.

Brightness mapping

The GUI version uses a slightly more advanced brightness mapping which provides more adjustability and better results. The GUI is set up to tweak settings and restack easily.

xp = [0, kneein, peak]
yp = [0, kneeout, 255]
img_array = np.interp(img_array, xp, yp).astype("uint8")

Colorize and stack

To stack the layers together, I use the colors from the GUI colorpicker to make a color image of each layer, then "stack" them with np.fmax.

# starting here with lists of img_arrays and color_rgb
img_stack = None
for i, img_array in enumerate(img_arrays):
    r = np.multiply(img_array, color_rgb[i][0] / 255.0)
    g = np.multiply(img_array, color_rgb[i][1] / 255.0)
    b = np.multiply(img_array, color_rgb[i][2] / 255.0)
    layer_color = np.dstack((r, g, b))
    if img_stack is None:
        img_stack = layer_color
    else:
        img_stack = np.fmax(layer_color, img_stack)
img_stack = img_stack.astype("uint8")

All of that seems to be working, but the results just don't compare:

In particular note how much larger the blown-out portion of the star in the center of mine is, and how the "textured" features of the nebula are much brighter in the NASA image.

Question 1: brightness mapping

I don't know what I don't know here. I think there's a lot of room for easy improvement through better brightness mapping, but it seems beyond just picking the right settings with my "single-knee" interpolation.

So my first question is are there other, better techniques for brightness mapping? In this dataset, the interesting nebula features which I'm not bringing out have pre-mapped brightness values in the 5-40 range, while the stars' brightness climbs quickly into the 100s and peaks in the multiple 1000s.

Adjusting the "knee" settings to get the star glare to match NASA and the nebula is barely visible, but adjust it so the nebula is as bright and the stars are way blown out. Would additional points in the interpolation help? If so, what's the strategy? Is there some topic or technique I just don't know about that I should google and read about that could help me develop here?

Question 2: layer stacking

My second question is am I even stacking the layers correctly? Basically I'm combining color images with np.fmax, taking the brightest pixel values of each channel from the stack. Because the layers are assigned arbitrary colors, not necessarily just primaries, summing them is definitely wrong but I'm wondering if fmax is also wrong. It certainly produces strange results, where layers can seem to obscure one another.

Question 3: other fits data

Last question, is there additional information in the fits file I should be using in some way, like to reduce the magnitude of the star brightness, or for anything else? There is a lot in there that I'm not using.