Stack multiple images in python pillow

Question

I am trying to do some image stacking in python pillow. What I would like to do is take a large number of images (say 10), and then for each pixel, take the median value like this: http://petapixel.com/2013/05/29/a-look-at-reducing-noise-in-photographs-using-median-blending/.

Right now, I can do it in an incredibly brute force manner (using getpixel and put pixel), but it takes a very long time.

Here is what I have so far:

import os
from PIL import Image
files = os.listdir("./")
new_im = Image.new('RGB', (4000,3000))
ims={}
for i in range(10,100):
    ims[i]={}           
    im=Image.open("./"+files[i])
    for x in range(400):
        ims[i][x]={}            
        for y in range(300):
            ims[i][x][y]=im.getpixel((x,y))

for x in range(400):
    for y in range(300):
        these1=[]
        these2=[]
        these3=[]
        for i in ims:
            these1.append(ims[i][x][y][0])
            these2.append(ims[i][x][y][1])
            these3.append(ims[i][x][y][2])
        these1.sort()
        these2.sort()
        these3.sort()
        new_im.putpixel((x,y),(these1[len(these1)/2],these2[len(these2)/2],these3[len(these3)/2]))

new_im.show()

Answer 1

You can vectorise a lot of these loops with arrays. For example np.array(im) will return an array of the pixels, with a shape (400, 300, 3). So to store everything in an array.

image_stacks = np.zeros(shape=(10, 400, 300, 3), dtype=np.uint8)
for i in xrange(image_stacks.shape[0]):
    # open image file and store in variable `im`, then
    image_stacks[i] = np.array(im)

Now you can calculate the median with your preferred way, but numpy has a method for that, too.

image_median = np.median(image_stacks, axis=0).astype(np.uint8)
image = Image.fromarray(image_median)

Two things to notice here is that np.median() will return a float type array, and we want to convert that to unsigned int8. The other thing is that if the number of elements is even, the median is calculated as the mean of the two middle values, which may end up being an odd number divided by two, such as 13.5. But when this is converted to integer, it will be rounded down. Such minor precision loss should not visually affect your result.