Travis Cunningham

Recently, I have been diving into astrophotography using a Nikon D3300 DSLR camera. Although the camera only has a 24.2 megapixel sensor, it performs well for beginners and I have gotten a lot of use out of it. Equatorial camera mounts are common, which are used to keep camera alignment on objects in the sky as the earth rotates. These are useful for long camera exposures which allow the camera to capture more light. However, equatorial mounts are not required to begin taking beautiful pictures of the night sky.

There are several options for processing astrophotography images if you don’t have an equatorial camera mount, such as Sequator or DeepSkyStacker. These software packages allow photographers to stack non-blurry star images with short exposures, which yields a result very similar to long-exposure images using an equatorial camera mount. Additionally, if star images are not aligned, these software tools will automatically align images and remove complex distortions caused by various camera lenses.

Today, I’ll be showing off my first iteration of a Jupyter notebook I’ve been using to process, align, and stack astrophotography images. This notebook is a work in progress and represents the introduction to a fully featured image processing pipeline.

To begin, we will take a look at the Orion Nebula (Messier 42 or M42). The Orion Nebula is a nebula in the constellation Orion, with an apparent visual magnitude of 4 and an angular size of 85×60 arc-minutes. Due to its size and brightness, it’s super easy to identify with the naked eye.

The goal is to build a pipeline to produce images similar to the following image. Here is a final image of the Orion Nebula I took and processed with Sequator with some color/contrast edits in Photoshop:

Starting Images

Using my Nikon D3300, I took roughly 100 photographs using a 300mm lens at 6400 ISO with 1 second exposures, adding up to almost 1 minute of exposure time. I also took 60 photographs using the same camera settings with the lens cap on, which we will use to subtract camera sensor noise from our light frames in a later step. This step will be crucial in eliminating noise that is present at high ISO settings.

Here is a single dark frame (noise only) and a single light frame (signal + noise) side by side:

The light frame is very noisy! Let’s clean it up.

Install Python Dependencies

First, here are the Python dependencies to install in your terminal:

pip3 install numpy pillow rawpy tqdm astroalign

Noise Removal

To remove the noise from my light frames, I combined all 60 dark frame images and saved the average pixel value for each pixel in the image. This forms an ‘average’ noise frame which will be used to subtract noise from the light frames. The dark frames are saved in the ‘darks’ directory:

import os

import numpy as np
from PIL import Image
import rawpy
from tqdm.notebook import tqdm


def open_numpy_raw(filename: str) -> np.array:
    """Open raw image files and convert to numpy array."""
    raw = rawpy.imread(filename)
    rgb = raw.postprocess()
    return rgb
   
   
def generate_average_dark(directory: str) -> np.array:
    """Generate an average dark image from a directory of dark images."""
    images = []
    
    # Load all images into np arrays
    pbar = tqdm(total=len(os.listdir(directory)), desc='Reading dark image files')
    for filename in os.listdir(directory):
        relative_path = os.path.join(directory, filename)
        np_image = open_numpy_raw(relative_path)
        images.append(np_image)
        pbar.update(1)
    pbar.close()

    # Generate a new image based on the average value of each pixel
    np_images = np.array(images)
    np_average = np.mean(np_images, axis=0)
    return np_average
    

average_dark_image = generate_average_dark('darks')

Next, I subtracted the pixel value of the average dark image from each light image. The light frames are saved in the ‘lights’ directory. I created the ‘tmp’ directory manually:

def remove_noise_from_light_image(light_image: np.array, dark_image: np.array) -> np.array:
    """Remove noise from a light image using a previously generated average dark image."""
    width = light_image.shape[0]
    height = light_image.shape[1]
    
    if dark_image.shape[0] != width or dark_image.shape[1] != height:
        raise Exception(f"Light and dark images have different dimensions: {dark_image.shape} and {light_image.shape}")
    
    new_image = np.subtract(light_image, dark_image)
    black_image = np.zeros([width, height, 3], dtype=np.uint8)
    return np.maximum(new_image, black_image)


def remove_noise_from_lights(directory: str, output_directory: str, dark_image: np.array) -> None:
    """Remove noise from all light images in a directory using a previously generated average dark image."""
    pbar = tqdm(total=len(os.listdir(directory)), desc=f'Removing noise from light images in directory: {directory}')
    for filename in os.listdir(directory):
        relative_path = os.path.join(directory, filename)
        light_image = open_numpy_raw(relative_path)
        denoised_image = remove_noise_from_light_image(light_image, dark_image)
        
        light_pil_image = Image.fromarray(denoised_image.astype("uint8"))
        light_pil_image.save(os.path.join(output_directory, f"denoised-{filename}.png"))
        
        pbar.update(1)
    pbar.close()


remove_noise_from_lights(directory='lights', output_directory='tmp', dark_image=average_dark_image)

The resulting image removes a lot of the noise introduced by the high ISO setting:

Image Stacking

There is still a lot of noise in this image, which is where stacking comes into play. Using the astroalign Python library, I aligned all of the images to be ready for stacking. Next, the average of each pixel was saved into a final image. All images in the ‘tmp’ directory are the targets to be processed. The ‘tmpstacked’ directory will be used for stacking aligned files.

from functools import lru_cache

import astroalign


@lru_cache(maxsize=1024)
def open_pil_img(path: str) -> Image:
    return Image.open(path)


@lru_cache
def listdir(path: str):
    return os.listdir(path)


def stack_and_align_light_images(directory: str, output_directory: str, star_detection_sigma: int) -> np.array:
    """Align and stack images in a source directory."""
    target_image_name = os.listdir(directory)[0]
    target_image_relative_path = os.path.join(directory, target_image_name)
    target_img = Image.open(target_image_relative_path)
    
    # Align images and save to output_directory
    pbar = tqdm(total=len(os.listdir(directory)), desc=f'Aligning images in directory: {directory}')
    for filename in os.listdir(directory):
        relative_path = os.path.join(directory, filename)
        source_img = Image.open(relative_path)
        
        # Align the image to the target image
        registered, footprint = astroalign.register(source_img, target_img, detection_sigma=star_detection_sigma)
        
        save_to = os.path.join(output_directory, filename)
        stacked_pil_image = Image.fromarray(registered.astype("uint8"))
        stacked_pil_image.save(save_to)
        
        pbar.update(1)
    pbar.close()
    
    # Stack images in output_directory
    width = target_img.width
    height = target_img.height
    new_image = np.zeros([width, height, 3], dtype=np.uint8)
    pbar = tqdm(total=width, desc=f'Stacking images in: {output_directory}')
    for x in range(width):
        for y in range(height):
            pixel_values = []
            for filename in listdir(output_directory):
                relative_path = os.path.join(output_directory, filename)
                source_img = open_pil_img(relative_path)
                pixel_values.append(source_img.getpixel((x, y)))
            pixel_values_np = np.array(pixel_values)
            avg_value = np.average(pixel_values_np, axis=0)
            new_image[x, y] = avg_value
        pbar.update(1)
    pbar.close()
    
    return new_image


stacked_image = stack_and_align_light_images(directory='tmp', output_directory='tmpstacker', star_detection_sigma=2)
# Convert to a PNG:
stacked_pil_image = Image.fromarray(stacked_image.astype("uint8"))
stacked_pil_image.save("stacked.png")

The resulting image is more clear and removes a lot of noise. Here is 60 frames stacked:

As you can see, the final image is more clear than what we started with, however, there is a lot of room for improvement.

If we stack 400 frames, the resulting image has even more clarity:

Gamma Adjustments

Next, we can apply gamma correction to bring out all colors in the image. As you can see from the image above, there is not much color to it. First, let’s create a histogram function to visualize the distribution of colors in the image:

%matplotlib inline

# pip3 install opencv-python opencv-contrib-python
# apt-get install ffmpeg libsm6 libxext6

import cv2
import numpy as np
from matplotlib import pyplot as plt

def generate_color_histogram(path: str) -> None:
    cv2_img = cv2.imread(path)
    colors = ('b','g','r')
    plt.figure()
    for i, color in enumerate(colors):
        histr = cv2.calcHist([cv2_img], [i], None, [256], [0, 256])
        plt.plot(histr ,color=color)
        plt.xlim([0,256])
    plt.show()

Next, we can apply gamma adjustments to increase the colors of the image.

%matplotlib inline
from matplotlib.pyplot import imshow


def gamma_adjust(gamma: int, path: str) -> np.array:
    np_image = imageio.imread(path)
    height = np_image.shape[0]
    width = np_image.shape[1]
    new_image = np.zeros([width, height, 3], dtype=np.uint8)
    
    pbar = tqdm(total=width, desc=f'Gamma correcting image: {path}')
    for x in range(width):
        for y in range(height):
            red, green, blue = np_image[y, x]
            new_red = pow(red / 255, (1 / gamma)) * 255
            new_green = pow(green / 255, (1 / gamma)) * 255
            new_blue = pow(blue / 255, (1 / gamma)) * 255
            new_color = (int(new_red), int(new_green), int(new_blue))
            new_image[x, y] = new_color
        pbar.update(1)
    pbar.close()
    return new_image


curve_stretched_image = gamma_adjust(gamma=2, path='stacked.png')
imshow(curve_stretched_image, aspect='auto')
curve_stretched_pil_image = Image.fromarray(curve_stretched_image.astype(np.uint8))
curve_stretched_pil_image.save("gamma-enhancement.png")

generate_color_histogram(path='gamma-enhancement.png')

The resulting image is much brighter, which enhances the colors in the nebula:

The resulting histogram is still quite narrow:

Color Levels Adjustments and Normalization

Next, we can ‘stretch’ and normalize the distribution of colors:

import imageio.v2 as imageio


def normalize_levels(path: str) -> np.array:
    np_image = imageio.imread(path)
    norm = (np_image - np.min(np_image)) / (np.max(np_image) - np.min(np_image))
    return norm * 255


normalized_levels_image = normalize_levels(path='gamma-enhancement.png')
imshow(normalized_levels_image / 255, aspect='auto')
normalized_pil_image = Image.fromarray(normalized_levels_image.astype(np.uint8))
normalized_pil_image.save("normalized-levels.png")

generate_color_histogram(path='normalized-levels.png')

The resulting image is much more colorful:

The resulting histogram shows a wider range of colors and a large amount of red and green:

To better balance the colors, we can adjust the colors by subtracting from red and green pixels. The r, g, and b parameters passed to the decrease_colors function were chosen by hand based on the previous histogram:

import imageio.v2 as imageio


def modify_colors(path: str, r: int, g: int, b: int) -> np.array:
    """Modify colors by a range from -255 to 255 on a per-color basis."""
    np_image = imageio.imread(path)
    width = np_image.shape[0]
    height = np_image.shape[1]
    new_image = np.zeros([width, height, 3], dtype=np.uint8)
    
    pbar = tqdm(total=width, desc=f'Color correcting image: {path}')
    for x in range(width):
        for y in range(height):
            red, green, blue = np_image[x, y]
            new_red = min(max(red + r, 0), 255)
            new_green = min(max(green + g, 0), 255)
            new_blue = min(max(blue + b, 0), 255)
            new_color = (int(new_red), int(new_green), int(new_blue))
            new_image[x, y] = new_color
        pbar.update(1)
    pbar.close()
    return new_image


corrected_colors_image = modify_colors(path='normalized-levels.png', r=-30, g=-30, b=0)
imshow(corrected_colors_image / 255, aspect='auto')
corrected_colors_pil_image = Image.fromarray(corrected_colors_image.astype(np.uint8))
corrected_colors_pil_image.save("corrected-colors.png")

generate_color_histogram(path='corrected-colors.png')

The resulting image has a better color balance:

The resulting histogram verifies the colors are more balanced:

There are several improvements I wish to make after this first iteration:

• Remove distortion effects caused by the camera lens, which is noticeable on the outside edges of the image.
• Perform curve-stretching, color levels, brightness, and saturation adjustments to bring out features of the nebula.
• Perform noise reduction and contrast adjustments
• Integrate with Starnet to generate starless images and edit star layers separately from the nebula

Stay tuned for more updates on my astrophotography adventure. If you have any comments or suggestions, please do so below.

Clear skys and happy hacking 🙂