Fast data structure or algorithm to find mean of each pixel in a stack of images

Question

I have a stack of images in which I want to calculate the mean of each pixel down the stack.

For example, let (x_n,y_n) be the (x,y) pixel in the nth image. Thus, the mean of pixel (x,y) for three images in the image stack is:

mean-of-(x,y) = (1/3) * ((x_1,y_1) + (x_2,y_2) + (x_3,y_3))

My first thought was to load all pixel intensities from each image into a data structure with a single linear buffer like so:

|All pixels from image 1| All pixels from image 2| All pixels from image 3|

To find the sum of a pixel down the image stack, I perform a series of nested for loops like so:

for(int col=0; col<img_cols; col++)
{
    for(int row=0; row<img_rows; row++)
    {
        for(int img=0; img<num_of_images; img++)
        {
            sum_of_px += px_buffer[(img*img_rows*img_cols)+col*img_rows+row];
        }
    }
}

Basically img*img_rows*img_cols gives the buffer element of the first pixel in the nth image and col*img_rows+row gives the (x,y) pixel that I want to find for each n image in the stack.

Is there a data structure or algorithm that will help me sum up pixel intensities down an image stack that is faster and more organized than my current implementation?

I am aiming for portability so I will not be using OpenCV and am using C++ on linux.

Answer 1

The problem with the nested loop in the question is that it's not very cache friendly. You go skipping through memory with a long stride, effectively rendering your data cache useless. You're going to spend a lot of time just accessing the memory.

If you can spare the memory, you can create an extra image-sized buffer to accumulate totals for each pixel as you walk through all the pixels in all the images in memory order. Then you do a single pass through the buffer for the division.

Your accumulation buffer may need to use a larger type than you use for individual pixel values, since it has to accumulate many of them. If your pixel values are, say, 8-bit integers, then your accumulation buffer might need 32-bit integers or floats.