How to do a median projection of a large image stack in Matlab

Question

I have a large stack of 800 16bit gray scale images with 2048x2048px. They are read from a single BigTIFF file and the whole stack barely fits into my RAM (8GB). Now I need do a median projection. That means I want to compute the median of each pixel across all 800 frames. The Matlab median function fails because there is not enough memory left make a copy of the whole array for the function call. What would be an efficient way to compute the median?

I have tried using a for loop to compute the median one pixel at a time, but this is still terribly slow.

Answer 1

Iterating over blocks, as @Shai suggests, may be the most straightforward solution. If you do have this problem frequently, you may want to consider converting the image to a mat-file, so that you can access the pixels as n-d array directly from disk.

%# convert to mat file
matObj = matfile('dest.mat','w');
matObj.data(2048,2048,numSlices) = 0; 
for t = 1:numSlices
    matObj.data(:,:,t) = imread(tiffFile,'index',t);
end

%# load a block of the matfile to take median (run as part of a loop)
medianOfBlock = median(matObj.data(1:128,1:128,:),3);

Answer 2

I bet that the distributions of the individual pixel values over the stack (i.e. the histograms of the pixel jets) are sparse.

If that's the case, the amount of memory needed to keep all the pixel histograms is much less than 2K x 2K x 64k: you can use a compact hash map to represent each histogram, and update them loading the images one at a time. When all updates are done, you go through your histograms and compute the median of each.