How to create a good image of "pillars of creation" by "make_lupton_rgb of astropy" in astropy,python?

Question

I am trying to write a code to generate an RBG image of the Pillars of Creation. For that I am using fits file corresponding to red, blue and green, and trying to use make_lupton_rbg to generate the RBG image. However I am getting full green image. I believe, I have to make adjustments to Q and stretch values, but I can't find anything to give it a good color (as seen in the pictures).

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.visualization import make_lupton_rgb

forc=np.float_()
r=fits.open("./673nmos.fits")[0].data
g=fits.open("./656nmos.fits")[0].data
b=fits.open("./502nmos.fits")[0].data

r = np.array(r,forc)
g = np.array(g,forc)
b = np.array(b,forc)

rgb_default = make_lupton_rgb(r,g,b,Q=1,stretch=0.1,filename="pillar.png")
plt.imshow(rgb_default, origin='lower')
plt.show()

The fits file were download from here

This is the output I am getting

And this is the output I should get (or at least something like it)

Answer 1

Scaling the r, g, and b arrays based on their relative brightnesses and using a high linear stretch gets much closer:

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.visualization import make_lupton_rgb

forc=np.float_()
r=fits.open("/path/to/673nmos/673nmos.fits")[0].data
g=fits.open("/path/to/656nmos/656nmos.fits")[0].data
b=fits.open("/path/to/502nmos/502nmos.fits")[0].data

r = np.array(r,forc)
g = np.array(g,forc)
b = np.array(b,forc)

rgb_default = make_lupton_rgb(r*5,g*0.75,b*8,Q=0.001,stretch=300,filename="pillar.png")
plt.imshow(rgb_default, origin='lower')
plt.show()

But clearly the linear stretching handles the spikes poorly, this can be compensated for by simply threshold filtering before applying make_lupton_rgb:

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.visualization import make_lupton_rgb

forc=np.float_()
r=fits.open("/path/to/673nmos/673nmos.fits")[0].data
g=fits.open("/path/to/656nmos/656nmos.fits")[0].data
b=fits.open("/path/to/502nmos/502nmos.fits")[0].data

r = np.array(r,forc)*5
g = np.array(g,forc)*0.75
b = np.array(b,forc)*8

t = 250
r[r > t] = t
g[g > t] = t
b[b > t] = t

rgb_default = make_lupton_rgb(r,g,b,Q=0.001,stretch=300,filename="pillar.png")
plt.figure(figsize=(8,8))
plt.imshow(rgb_default, origin='lower')
plt.show()