Stacking astronomy images with Python

Question

I thought this was going to be easier but after a while I'm finally giving up on this, at least for a couple of hours...

I wanted to reproduce this a trailing stars image from a timelapse set of pictures. Inspired by this:

The original author used low resolution video frames taken with VirtualDub and combined with imageJ. I imagined I could easily reproduce this process but with a more memory-conscious approach with Python, so I could use the original high-resolution images for a better output.

My algorithm's idea is simple, merging two images at a time, and then iterating by merging the resulting image with the next image. This done some hundreds of times and properly weighing it so that every image has the same contribution to the final result.

I'm fairly new to python (and I'm no professional programmer, that'll be evident), but looking around it appears to me the Python Imaging Library is very standard, so I decided to use it (correct me if you think something else would be better).

Here's what I have so far:

#program to blend many images into one
import os,Image
files = os.listdir("./")
finalimage=Image.open("./"+files[0]) #add the first image
for i in range(1,len(files)): #note that this will skip files[0] but go all the way to the last file
  currentimage=Image.open("./"+files[i])
  finalimage=Image.blend(finalimage,currentimage,1/float(i+1))#alpha is 1/i+1 so when the image is a combination of i images any adition only contributes 1/i+1.
  print "\r" + str(i+1) + "/" + str(len(files)) #lousy progress indicator
finalimage.save("allblended.jpg","JPEG")

This does what it's supposed to but the resulting image is dark, and if I simply try to enhance it, it's evident that information was lost due lack of depth in pixel's values. (I'm not sure what the proper term here is, color depth, color precision, pixel size). Here's the final result using low resolution images:

or one I was trying with the full 4k by 2k resolution (from another set of photos):

So, I tried to fix it by setting the image mode:

firstimage=Image.open("./"+files[0])
size = firstimage.size
finalimage=Image.new("I",size)

but apparently Image.blend does not accept that image mode.

ValueError: image has wrong mode

Any ideas?

(I also tried making the images "less dark" by multiplying it before combining them with im.point(lambda i: i * 2) but results were just as bad)

Answer 1

The problem here is that you are averaging the brightness at each pixel. This may seem sensible but it is actually not what you want at all -- the bright stars will get "averaged away" because they move accross the image. Take the following four frames:

1000 0000 0000 0000
0000 0100 0000 0000
0000 0000 0010 0000
0000 0000 0000 0001

If you average those, you will get:

0.25 0    0    0
0    0.25 0    0
0    0    0.25 0
0    0    0    0.25

When you want:

1000
0100
0010
0001

Instead of blending the images you can try taking the maximum seen in any image for each pixel. If you have PIL you can try the lighter function in ImageChops.

from PIL import ImageChops
import os, Image
files = os.listdir("./")
finalimage=Image.open("./"+files[0])
for i in range(1,len(files)):
    currentimage=Image.open("./"+files[i])
    finalimage=ImageChops.lighter(finalimage, currentimage)
finalimage.save("allblended.jpg","JPEG")

Here is what I got:

EDIT: I read the Reddit post and see that he actually combined two approaches -- one for the star trails and a different one for the Earth. Here is a better implementation of the averaging you tried, with proper weighting. I used a numpy array for the intermediate storage instead of the uint8 Image array.

import os, Image
import numpy as np
files = os.listdir("./")
image=Image.open("./"+files[0])
im=np.array(image,dtype=np.float32)
for i in range(1,len(files)):
    currentimage=Image.open("./"+files[i])
    im += np.array(currentimage, dtype=np.float32)
im /= len(files) * 0.25 # lowered brightness, with magic factor
# clip, convert back to uint8:
final_image = Image.fromarray(np.uint8(im.clip(0,255)))
final_image.save('all_averaged.jpg', 'JPEG')

Here is the image, which you could then combine with the star trails from the previous one.