image stack population is slow in numpy
Question: I am reading stack of separate tiff's into single 3D array via numpy/python.
When files are just read and plugged into some variable, speed scales linearly
with number of files, for example, loading 100 files takes 0.2s, loading 1000
files takes 2.46s and so on.

However, when I try to create a 3D stack out of those files, using _dstack()_
time starts to scale non-linearly, e.g. 0.21 sec for 10 files, 5.39sec for 100
files and so on.

I realize that slow-down is caused by some magic behind dstack(). What is the
proper and fastest way to create 3D stacks from set of image files?

If I don't use dtack(), but pre-create 3D array an then populate it, script is
running faster, but still scales non-linearly (0.2s for 10, 2.2s for 100, 40s
for 1000 images) Code for this case:

    
    
    import numpy as np
    from PIL import Image
    import time
    import random
    
    def toc(t):
        return time.time() - t
    
    i_max = 1000
    t = time.time()
    for i in range(0,i_max):
        fname = r"..\Pos0\img_"+("%09d"%i)+"_Default_000.tif"
        im = np.array(Image.open(fname))
        if i>0:
            stack[:,:,i] = im
        else:
            s = im.shape
            stack = np.empty((s[0],s[1],i_max))
            stack[:,:,0] = im
    print toc(t)
    

PS: Python 2.7.8 Anaconda 2.1, Intel i5 @ 32GB RAM, reading from 4-striped
HDDs


Answer: You can solve this problem by changing the order of your indexing, making the
image index first. Like this:

    
    
    i_max = 1000
    sx, sy = 1000,1000
    t = time.time()
    for i in range(0,i_max):
        im = np.ones((sx,sy))
        if i>0:
            #stack[:,:,i] = im
            stack[i,:,:] = im
        else:
            #stack = np.empty((sx,sy,i_max))
            #stack[:,:,0] = im
            stack = np.empty((i_max, sx, sy))
            stack[0,:,:] = im
    print toc(t)
    

The times I get are:

    
    
    4.44851183891   # original order, i_max = 100
    118.510767937   # original order, i_max = 1000
    
    1.78239989281   # modified order, i_max = 100
    23.4904351234   # modified order, i_max = 1000
    

I think the reason for this is that the data is not handled bit-by-bit, but in
chunks of fixed size, regardless of how much data in each chunk is actually
being used. With your original order, the chunks were not being used
efficiently. For example, see the first few minutes of this video:
<https://vimeo.com/97337258> That is, it's not about RAM but CPU cache.

