Converting to RGB

If you have two or more images of the same scene taken at different wavelengths, you may wish to combine them to create a colour composite.

For ultimate control, you can do this manually using imview. Simply map your channels to RGB values using imview and then sum the results.

For convenience, AstroImages.jl provides the function composecolors.

Using composecolors

We'll demonstrate composecolors using Hubble images of the Antenae colliding galaxies.

One can be very scientific about this process, but often the goal of producing color composites is aesthetic or about highlighting certain features for public consumption.

I'll set the default color map to grayscale to avoid confusion.

using AstroImages

AstroImages.set_cmap!(nothing)

Let's start by downloading the separate color channel FITS files:

antred = AstroImage(download("http://www.astro.uvic.ca/~wthompson/astroimages/fits/antenae/red.fits"))

antgreen = AstroImage(download("http://www.astro.uvic.ca/~wthompson/astroimages/fits/antenae/green.fits"))

antblue = AstroImage(download("http://www.astro.uvic.ca/~wthompson/astroimages/fits/antenae/blue.fits"))

anthalph = AstroImage(download("http://www.astro.uvic.ca/~wthompson/astroimages/fits/antenae/hydrogen.fits")); # Hydrogen-Alpha; we'll revisit later

The images will have to be aligned and cropped to the same size before making a color composite.

In order to compose these images, we'll have to match the relative intensity scales and clip outlying values. Thankfully, composecolors handles most of these details automatically.

rgb1 = composecolors([antred, antgreen, antblue])

It's a start!

By default, if you provide three images these are mapped to the color channels red, green, and blue. The intensities are limited to Percent(99.5).

We can now tweak these defaults to our tastes. We could try clamping the intensities more agressively to bring out more of the galaxy structure:

rgb2 = composecolors(
    [antred, antgreen, antblue],
    clims=Percent(97)
)

This looks okay but saturates the galaxy cores.

Let's take care of that gash through the image by just blanking it out.

mask = antgreen .== antgreen[end,begin]
# remove holes in the mask
using ImageFiltering, Statistics
mask = BitMatrix(mapwindow(median, mask, (3,3)))
imview(mask)

antred[mask] .= NaN
antgreen[mask] .= NaN
antblue[mask] .= NaN
anthalph[mask] .= NaN

Typically we need to perform a "gamma correction" aka non-lienar stretch to map the wide dynamic range of astronomical images into a narrower human visible range. We can do this using the stretch keyword. An asinhstretch is typically recommended when preparing RGB images:

rgb3 = composecolors(
    [antred, antgreen, antblue],
    stretch=asinhstretch
)

Keywords like strech, clims, etc can be either a single value for all channels or a list of separate values/functions per channel.

The green channel appears to be quite faint compared to the red and blue channels. We can modify that by adjusting the relative intensities of the channels.

We could also do this using a combination of the contrast and bias keywords.

rgb4 = composecolors(
    [antred, antgreen, antblue],
    stretch=asinhstretch,
    multiplier=[1,1.7,1]
)

That's better! Let's go one step further, and incorporate a fourth chanel: Hydrogen Alpha. Hydrogen Alpha is a narrow filter centered around one of the emission lines of Hydrogen atoms. It traces locations with hot gas; mostly star-formation regions in this case.

imview(anthalph, cmap=:magma, clims=Zscale())

We'll now need to specify the color channels we want to use for each wavelength since we can't use just the default three RGB. We can use any named color or julia ColorScheme.

rgb5 = composecolors(
    [antred, antgreen, antblue, anthalph],
    ["red", "green", "blue", "maroon1"],
    stretch=asinhstretch,
    multiplier=[1,1.7,1,0.8]
)

Additionally, I'd like to just show the brightest areas of Hydrogen alpha emission rather than adding a diffuse pink glow. We can turn off the stretch for this one channel:

rgb6 = composecolors(
    [antred, antgreen, antblue, anthalph],
    ["red", "green", "blue", "maroon1"],
    stretch=[
        asinhstretch,
        asinhstretch,
        asinhstretch,
        identity,
    ],
    multiplier=[1,1.7,1,0.8]
)

Finally, we can crop the image and save it as a PNG.

crop = rgb6[200:end-100,50:end-50]

save("antenae-composite.png", crop)

If you want to save it in a format like JPG that doesn't support transparent pixels, you could replace the masked area with zeros instead of NaN.