extracting many regions of interests ROIs) from thousand images

Question

I have a large set of microscopy images and each image has several hundreds of spots (ROIs). These spots are fixed in space. I want to extract each spot from each image and save into workspace so that I can analyze them further.

I have written a code myself and it working perfectly but its too slow. It takes around 250 sec to completely read out all the spots from every image.

The core of my code looks as following:

for s=1:NumberImages   
  im1=imread(fn(s,1).name);    
  im=im1-medfilt2(im1,[15,15]);    
  for i=1:length(p_g_x)    
    GreenROI(i,s)=double(sum(sum(im(round(p_g_y(i))+(-2:2), round(p_g_x(i))+(-2:2)))));
    RedROI(i,s)=double(sum(sum(im(round(p_r_y(i))+(-2:2), round(p_r_x(i))+(-2:2)))));        
  end
end

As you can see from the code I am extracting 5x5 regions. Length of p_g_x is between 500-700.

Thanks for your input. I used profile viewer to figure out which function exactly is taking more time. It was median filter which is taking a lot of time (~90%).

Any suggestion to fast it up will be greatly appreciated.

thanks

Mahipal

Answer 1

Vectorized code that operates on each image -

%// Pre-compute green and red indices to be used across all the images
r1 = round(bsxfun(@plus,permute(p_g_y,[3 2 1]),[-2:2]'));
c1 = round(bsxfun(@plus,permute(p_g_x,[3 2 1]),[-2:2]));
green_ind = reshape(bsxfun(@plus,(c1-1)*size(im,1),r1),[],numel(p_g_x));

r2 = round(bsxfun(@plus,permute(p_r_y,[3 2 1]),[-2:2]'));
c2 = round(bsxfun(@plus,permute(p_r_x,[3 2 1]),[-2:2]));
red_ind = reshape(bsxfun(@plus,(c2-1)*size(im,1),r2),[],numel(p_g_x));

for s=1:NumberImages
    im1=imread(fn(s,1).name);
    im=double(im1-medfilt2(im1,[15,15]));

    GreenROI=sum(im(green_ind));
    RedROI =sum(im(red_ind));
end

Answer 2

Use Matlab's profiling tools!

profile on % Starts the profiler
% Run some code now.
profile viewer % Shows you how often each function was called, and
               % where most time was spent. Try to start with the slowest part.
profile off  % Resets the Profiler, so you can measure again.

Pre-allocate

Preallocate the output because you know the size and this way it is much faster. (Matlab told you this already!)

GreenROI = zeros(length(p_g_x), NumberImages); % And the same for RedROI.

Use convolution

Read about Matlab's conv2 code.

for s=1:NumberImages   
  im=imread(fn(s,1).name);    
  im=im-medfilt2(im,[15,15]);    
  % Pre-compute the sums first. This will only be faster for large p_g_x
  roi_image = conv2(im, ones(5,5));
  for i=1:length(p_g_x)    
    GreenROI(i,s)=roi_image(round(p_g_y(i)), round(p_g_x(i))); % You might have to offset the indices by 2, because of the convolution. Check that results are the same.
    RedROI(i,s)=roi_image(round(p_r_y(i)), round(p_r_x(i)));        
  end
end

Matlab-ize the code

Now, that you've used convolution to get an image of sums over 5x5 windows (or you could've used @Shai's accumarray, same thing), you can speed things up further by not iterating through each element in p_g_x but use it as a vector straight away.

I leave that as an exercise for the reader. (convert p_g_x and p_g_y to indices using sub2ind, as a hint).

Update

Our answers, mine included, showed how premature optimisation is a bad thing. Without knowing, I assumed that your loop would take most of the time, but when you measured it (thanks!) it turns out that is not the problem. The bottleneck is medfilt2 the median filter, which takes 90% of the time. So you should address this first. (Note, that on my computer your original code is fast enough for my taste but it is still the median filter taking up most of the time.)

Looking at what the median filter operation does might help us figure out how to make it faster. Here is an image. On the left you see the original image. In the middle the median filter and on the right there is the result.

To me the result looks awfully similar to an edge detection result. (Mathematically this is no surprise.)

I would suggest you start experimenting with various edge detections. Have a look at Canny and Sobel. Or just use conv2(image, kernel_x) where kernel_x = [1, 2, 1; 0, 0, 0; -1, -2, -1] and the same but transposed for a kernel_y. You can find various edge detection options here: edge(im, option). I tried all options from {'sobel', 'canny', 'roberts', 'prewitt'}. Except for Canny, they all take about the same time as your median filter method. Canny is 4x slower, the rest (including the original) take 7.x seconds. All of this without a GPU. imgradient was 9 seconds.

So from this I would say that you can't get any faster. If you have a GPU and it works with Matlab, you could speed it up. Load your image as gpuArrays. There is an example on the medfilt2 documentation. You can still do minor speed ups but they can only amount to a 10% speed increase, so are hardly worthwile.

Answer 3

A few things you should do

1. Pre-allocate as suggested by Didac Perez.

2. Use profiler to see what exactly takes long in your code, is it the median filter? is it the indexing?

3. Assuming all images are of the same size, you can use accumarray and a fixed mask subs to quickly sum the values:

   subs_g = zeros( h, w ); %// allocate mask for green
   subs_r = zeros( h, w );  
   subs_g( sub2ind( [h w], round(p_g_y), round(p_g_x) ) = 1:numel(p_g_x); %//index each region
   subs_g = conv2( subs_g, ones(5), 'same' );
   subs_r( sub2ind( [h w], round(p_r_y), round(p_r_x) ) = 1:numel(p_r_x); %//index each region
   subs_r = conv2( subs_r, ones(5), 'same' );
   sel_g = subs_g > 0;
   sel_r = subs_r > 0;
   subs_g = subs_g(sel_g);
   subs_r = subs_r(sel_r);
   

   once these masks are fixed, you can process all images

   %// pre-allocation goes here - I'll leave it to you
   for s=1:NumberImages   
       im1=imread(fn(s,1).name);    
       im=double( im1-medfilt2(im1,[15,15]) ); 
       accumarray( subs_g, im( sel_g ) ); % summing all the green ROIs 
       accumarray( subs_r, im( sel_r ) ); % summing all the green ROIs 
   end
   

Answer 4

First, preallocate your GreenROI and RedROI structures since you previously know the final size. Now, you are resizing them again and again in each iteration.

Secondly, I do recommend you to use "tic" and "toc" to investigate where is the problem, it will give you useful timings.