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 š
Leave a Reply