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.
