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: ![Inspiration](http://i.stack.imgur.com/tZ7zW.jpg)

[The original
author](http://www.reddit.com/r/pics/comments/pkycp/star_trails_above_earthly_light_trails_as_viewed/c3qa141)
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](http://eol.jsc.nasa.gov/Videos/CrewEarthObservationsVideos/) 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:

![Low resolution result](http://i.stack.imgur.com/ixppn.jpg)

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

![High resolution result with another set of
images](http://i.stack.imgur.com/JbNWL.jpg)

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: 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: ![Low res image set
stacked](http://i.stack.imgur.com/esgaC.jpg)

**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. ![Low res images added
together](http://i.stack.imgur.com/RxrSU.jpg)

