python - Fastest approach to read thousands of images into one big numpy array

Question

I'm trying to find the fastest approach to read a bunch of images from a directory into a numpy array. My end goal is to compute statistics such as the max, min, and nth percentile of the pixels from all these images. This is straightforward and fast when the pixels from all the images are in one big numpy array, since I can use the inbuilt array methods such as .max and .min, and the np.percentile function.

Below are a few example timings with 25 tiff-images (512x512 pixels). These benchmarks are from using %%timit in a jupyter-notebook. The differences are too small to have any practical implications for just 25 images, but I am intending to read thousands of images in the future.

# Imports
import os
import skimage.io as io
import numpy as np

1. Appending to a list

   %%timeit
   imgs = []    
   img_path = '/path/to/imgs/'
   for img in os.listdir(img_path):    
       imgs.append(io.imread(os.path.join(img_path, img)))    
   ## 32.2 ms ± 355 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
   

2. Using a dictionary

   %%timeit    
   imgs = {}    
   img_path = '/path/to/imgs/'    
   for img in os.listdir(img_path):    
       imgs[num] = io.imread(os.path.join(img_path, img))    
   ## 33.3 ms ± 402 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
   

For the list and dictionary approaches above, I tried replacing the loop with a the respective comprehension with similar results time-wise. I also tried preallocating the dictionary keys with no significant difference in the time taken. To get the images from a list to a big array, I would use np.concatenate(imgs), which only takes ~1 ms.

3. Preallocating a numpy array along the first dimension

   %%timeit    
   imgs = np.ndarray((512*25,512), dtype='uint16')    
   img_path = '/path/to/imgs/'    
   for num, img in enumerate(os.listdir(img_path)):    
       imgs[num*512:(num+1)*512, :] = io.imread(os.path.join(img_path, img))    
   ## 33.5 ms ± 804 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
   

4. Preallocating a numpy along the third dimension

   %%timeit    
   imgs = np.ndarray((512,512,25), dtype='uint16')    
   img_path = '/path/to/imgs/'    
   for num, img in enumerate(os.listdir(img_path)):    
       imgs[:, :, num] = io.imread(os.path.join(img_path, img))    
   ## 71.2 ms ± 2.22 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
   

I initially thought the numpy preallocation approaches would be faster, since there is no dynamic variable expansion in the loop, but this does not seem to be the case. The approach that I find the most intuitive is the last one, where each image occupies a separate dimensions along the third axis of the array, but this is also the slowest. The additional time taken is not due to the preallocation itself, which only takes ~ 1 ms.

I have three question regarding this:

1. Why is the numpy preallocation approaches not faster than the dictionary and list solutions?
2. Which is the fastest way to read in thousands of images into one big numpy array?
3. Could I benefit from looking outside numpy and scikit-image, for an even faster module for reading in images? I tried plt.imread(), but the scikit-image.io module is faster.

Answer 1

Part A : Accessing and assigning NumPy arrays

Going by the way elements are stored in row-major order for NumPy arrays, you are doing the right thing when storing those elements along the last axis per iteration. These would occupy contiguous memory locations and as such would be the most efficient for accessing and assigning values into. Thus initializations like np.ndarray((512*25,512), dtype='uint16') or np.ndarray((25,512,512), dtype='uint16') would work the best as also mentioned in the comments.

After compiling those as funcs for testing on timings and feeding in random arrays instead of images -

N = 512
n = 25
a = np.random.randint(0,255,(N,N))

def app1():
    imgs = np.empty((N,N,n), dtype='uint16')
    for i in range(n):
        imgs[:,:,i] = a
        # Storing along the first two axes
    return imgs

def app2():
    imgs = np.empty((N*n,N), dtype='uint16')
    for num in range(n):    
        imgs[num*N:(num+1)*N, :] = a
        # Storing along the last axis
    return imgs

def app3():
    imgs = np.empty((n,N,N), dtype='uint16')
    for num in range(n):    
        imgs[num,:,:] = a
        # Storing along the last two axes
    return imgs

def app4():
    imgs = np.empty((N,n,N), dtype='uint16')
    for num in range(n):    
        imgs[:,num,:] = a
        # Storing along the first and last axes
    return imgs

Timings -

In [45]: %timeit app1()
    ...: %timeit app2()
    ...: %timeit app3()
    ...: %timeit app4()
    ...: 
10 loops, best of 3: 28.2 ms per loop
100 loops, best of 3: 2.04 ms per loop
100 loops, best of 3: 2.02 ms per loop
100 loops, best of 3: 2.36 ms per loop

Those timings confirm the performance theory proposed at the start, though I expected the timings for the last setup to have timings in between the ones for app3 and app1, but maybe the effect of going from last to the first axis for accessing and assigning isn't linear. More investigations on this one could be interesting (follow up question here).

To claify schematically, consider that we are storing image arrays, denoted by x (image 1) and o (image 2), we would have :

App1 :

[[[x 0]
  [x 0]
  [x 0]
  [x 0]
  [x 0]]

 [[x 0]
  [x 0]
  [x 0]
  [x 0]
  [x 0]]

 [[x 0]
  [x 0]
  [x 0]
  [x 0]
  [x 0]]]

Thus, in memory space, it would be : [x,o,x,o,x,o..] following row-major order.

App2 :

[[x x x x x]
 [x x x x x]
 [x x x x x]
 [o o o o o]
 [o o o o o]
 [o o o o o]]

Thus, in memory space, it would be : [x,x,x,x,x,x...o,o,o,o,o..].

App3 :

[[[x x x x x]
  [x x x x x]
  [x x x x x]]

 [[o o o o o]
  [o o o o o]
  [o o o o o]]]

Thus, in memory space, it would be same as previous one.

---

Part B : Reading image from disk as arrays

Now, the part on reading image, I have seen OpenCV's imread to be much faster.

As a test, I downloaded Mona Lisa's image from wiki page and tested performance on image reading -

import cv2 # OpenCV

In [521]: %timeit io.imread('monalisa.jpg')
100 loops, best of 3: 3.24 ms per loop

In [522]: %timeit cv2.imread('monalisa.jpg')
100 loops, best of 3: 2.54 ms per loop