Fast data structure or algorithm for finding the average value of each pixel in the image stack

Question

I have a stack of images in which I want to calculate the average value of each pixel on the stack.

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

 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 the pixel intensities from each image into a data structure using one linear buffer as follows:

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

To find the sum of a pixel in the image stack, I do a series of nested loops like:

 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 pixel (x,y) that I want to find for each n-image on the stack.

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

I am focused on portability, so I will not use OpenCV and use C ++ for Linux.

Answer 1

The problem with the nested loop in the question is that it is not very convenient for the cache. You will skip the memory skip in large steps, effectively making your data cache useless. You will spend a lot of time accessing memory.

If you can save memory, you can create an additional image size buffer to accumulate totals for each pixel when you go through all the pixels in all the images in memory. Then you do one pass through the division buffer.

You may need a larger type in your cumulative buffer than for individual pixel values, since it must accumulate many of them. If your pixel values ​​are, say, 8-bit integers, then your accumulation buffer may need 32-bit integers or floating elements.

Answer 2

Usually a stack of pixels

 (x_1,y_1),...,(x_n,y_n) 

conditionally independent of the stack

 (a_1,b_1),...,(a_n,b_n) 

And even if they were not (assuming a certain set of data), then modeling their interactions is a difficult task and will give you only an estimate of the average. So, if you want to calculate the exact average for each stack, you have no choice but to iterate through the three loops you supply. Languages ​​like Matlab / octave and libraries like Theano (python) or Torch7 (lua) parallelize these iterations. If you use C ++, then what you do is well suited for Cuda or OpenMP. As for portability, I think OpenMP is an easier solution.

Answer 3

A portable, fast data structure specifically for the average calculation can be:

 std::vector<std::vector<std::vector<sometype> > > VoVoV; VoVoV.resize(img_cols); int i,j; for (i=0 ; i<img_cols ; ++i) { VoVoV[i].resize(img_rows); for (j=0 ; j<img_rows ; ++j) { VoVoV[i][j].resize(num_of_images); // The values of all images at this pixel are stored continguously, // therefore should be fast to access. } } VoVoV[col][row][img] = foo; 

As a side note, 1/3 in your example will evaluate to 0 , which is not what you want.

To quickly add / average, you can now:

 sometype sum = 0; std::vector<sometype>::iterator it = VoVoV[col][row].begin(); std::vector<sometype>::iterator it_end = VoVoV[col][row].end(); for ( ; it != it_end ; ++it) sum += *it; sometype avg = sum / num_of_images; // or similar for integers; check for num_of_images==0 

Basically, you should not rely on the fact that the compiler will optimize the recalculation of always the same offsets.