any ideas how to make the code more efficient? Now it takes a few days!

I have 200 images and I am trying to calculate Signal to Noise ratio at each pixel following the formula

SNR at each pixel = mean(pixel (x,y) value from 200 image)/std(pixel (x,y) value from 200 image).

It takes ages to process the images using the code bellow even on super powerful computer. Any ideas as to how to speed up the code are very much appreciated.

% Reading multiple image files test1_0001.edf to test1_0200.edf in a loop

for i=1:200
    if i<100
        filename=strcat('test1_00',num2str(i),'.edf');
        [header matrix]=pmedf_read(filename);
        variablename=strcat('test1_00',num2str(i));
        eval([variablename '=matrix;']);
    else
       filename=strcat('test1_0',num2str(i),'.edf');
        [header matrix]=pmedf_read(filename);
        variablename=strcat('test1_0',num2str(i));
        eval([variablename '=matrix;']);
    end
end

% Getting mean count per pixel from all images
x=2500;
y=1500;
for r=1:x % rows
for c=1:y % columns
for i=1:200
if i<100
pixel_image(i)=eval(strcat('test1_00',num2str(i),'(r,c)'));
else
signal_value(i)=eval(strcat('test1_0',num2str(i),'(r,c)'));
end
signal_pixel(r,c)=sum(signal_value(:))/200;
std_pixel(r,c)=std(signal_value(:));
snr_pixel(r,c)=signal_pixel(r,c)/std_pixel(r,c);

end
end
end

[Cedric]

Di you profile the code to find the bottleneck(s)? You are doing several operations that can be slow (reading files, using EVAL, etc), and it would be relevant to know how they compare in term of execution time.

If you didn't profile it yet, modify your script so it treats e.g. 5 files instead of 200, then type the following in the command line:

 profile viewer

In the field, type the name of your script and click on [Start profiling]. Then study the report.

[Alex]

when I was trying to process 200 images the computer died and I had to reboot the station. Now I am trying to process just 5 images to see how long it takes. I started more than 7 hr ago and it is still running... As for the reading the images it is quite fast. But the calculations at each pixel take ages.

Answer 1

[Cedric]

You probably want to do something like the following (to adapt to your case, and check that it is working):

 n_img = 200 ;
 dim_x = 2500 ;     
 dim_y = 1500 ;

 buffer = zeros(dim_x, dim_y, n_img) ;

 for k = 1 : n_img
    filename = sprintf('test1_%03d.edf', k) ;
    [header, matrix] = pmedf_read(filename) ;
    buffer(:,:,k) = matrix ;
 end

 signal = mean(buffer, 3) ;
 noise  = std(buffer, [], 3) ;
 snr    = signal ./ noise ;

What we do here is to store each image as a "page" in a 3D array. We then operate along the 3rd dimension to compute the mean and std.

If you want to experiment to understand the mechanism on a simpler case, you can proceed as follows: generate two fake 3x4 "images" with random pixel values and compute mean, std, etc along the 3rd dim:

 >> buffer = randi(10, [3,4,2])
 buffer(:,:,1) =                           % Page/image 1.
     8     8    10     9
     7     7     7     4
     7    10     9     2
 buffer(:,:,2) =                           % Page/image 2
     8     7     2     2
     3     2     3     9
     6     3     9     3

 >> mean(buffer, 3)
 ans =
    8.0000    7.5000    6.0000    5.5000
    5.0000    4.5000    5.0000    6.5000
    6.5000    6.5000    9.0000    2.5000

 >> std(buffer, [], 3)
 ans =
         0    0.7071    5.6569    4.9497
    2.8284    3.5355    2.8284    3.5355
    0.7071    4.9497         0    0.7071

[Alex]

it is a very good idea to work with images in 3D. But I need SNR at each pixel (not just for the whole image) and it makes everything very complex. If we have 200 images then Signal(x1,y1) = ((1st image pixel value(x1,y1)+...200th pixel image value(x1/y1))/200 and Noise(x1,y1) = std(pixel value(x1,y1) across 200 image). Thus SNR at each pixel = mean(x,y)/std(x,y).

[Cedric]

This is precisely what we are doing here. Experiment with the example that I give at the end of my post. We operate along the 3rd dimension means that e.g. signal(x,y) is defined by the mean of pixel (x,y) over all images. If you look at the example, element (2,1) of the mean (=5) is the mean of element (2,1) of page/image 1 (=7) and element (2,1) of page/image 2 (=3). So you have to imaging e.g. that we stack flat images vertically, and that we perform computations along the vertical axis.

[Alex]

sorry, got confused. You are absolutely right and I am gonna test it right away.

Answer 2

[Walter Roberson]

The eval() are slowing you down.

http://matlab.wikia.com/wiki/FAQ#How_can_I_create_variables_A1.2C_A2.2C....2CA10_in_a_loop.3F