How to rewrite a loop in a logical way?

I have this loop but it takes too much time for processing it because of the big amount of elements of "x" and "y"... I have tried to write it in a more compact way but it was not successful... I have an image "im_raw" with "2546x3318" elements. It has a lot of holes which are set with "nan" values and I need to fill these "nan" values with the maximum of a submatrix around each gridpoint (9x9 neighbours). So where I have stopped is in this point when I have to find for each gridpoint, which must coincide with the x and y which contain the coordinates where the "nan" values are located and I want to replace by the maximum of that submatrix (9x9) centered in each gridpoint. Could somebody take a look at this example and send me a feedback? I would appreciate it... Thanks in advance.

for i = 10:max(x(:,1))-9

    for j = 10:max(y(:,1))-9

        im_raw(x(i,1),y(j,1)) = max(max(im_raw(i-9:1:i+9,j-9:1:j+9))) 

    end

end

[Stephen23]

Note that you should not use i or j as variable names, as these are the names of the inbuilt imaginary unit .

Answer 1

[Guillaume]

Is your neighbourhood 9x9 as in your description or 19x19 as in your code?

I'm sure there must be something for this in the image processing toolbox but I can't think of it right now. In any case, does this work for you?

im_raw(isnan(im_raw)) = -inf; %replace nans with -Inf so it doesn't mess up |max|
%pad matrix to avoid edge issues:
im_raw2 = ones(size(im_raw) + 18) * -Inf; %assuming a 19x19 neighbourhood
im_raw2(10:end-9, 10:end-9) = im_raw;   %assuming a 19x19 neighbourhood

[r, c] = find(isinf(im_raw)); %find values to replace
for p = 1:numel(r)
    neighbourhood = im_raw2(r(p):r(p)+18, c(p):c(p)+18);
    im_raw(r(p), c(p)) = max(neighbourhood(:));
end

edit: fixed bug pointed out by Michael

[Antonio]

Yes )x) I mean on each side of the central pixel, therefore is 19x19 or a radius of 9 pixels... Noe it is running but I knew it would take a lot of time, too much time... I will try your option and I will tell you... thanks for the fast feedback...

[Antonio]

Hello, after checking the result I have to tell you that first the image was too big and after the weekend processing it was still on processing... I stopped it and I built other image smaller just to see the effect of the code... I have seen that the "im_raw" matrix take only the last maximum of the last neighbourhood and build the image with theis maximum value for all the "nan" or "Inf" available values... Of course I would like to see that for every neighbourhood the corresponding maximum is stored in the final "im_raw" matrix and continues for the other "Inf" values to replace them with their corresponding maximum values of their neighbourhoods... Anyway, thanks...

[Michael Haderlein]

Guillaume's code has a little bug, I think. The line

    neighbourhood = im_raw2(r:r+19, c:c+19);

should rather be

    neighbourhood = im_raw2(r(p):r(p)+19, c(p):c(p)+19);

This should fix the issue you report.

Adding a -inf (or nan) frame around the image to avoid edge issues is perfectly fine, it just needs quite some extra memory as it copies your image. In case you are limited by memory, you could also keep the unframed image and use as neighbourhood:

    neighbourhood=img_raw(max(1,r(p)-9):min(end,r(p)+9),max(1,c(p)-9):min(end,c(p)+9));

A final comment on the very first line: Instead of overwriting all nans with -inf, you can also keep the nans and use nanmax instead of max. Though I have to admit that I didn't benchmark the two possibilities to check which one is faster.

[Guillaume]

Michael, thanks for the correction. Indeed, I forgot the indexing. I did test the code, it must be a line I didn't copy/paste. edit: and actually, the bounds are also wrong.

nanmax is indeed another option, assuming you have the stats toolbox.

Antonio, sorry about the bugs. I did test the code on a random image of the same size as yours, and it was fairly fast.

[Michael Haderlein]

Oh, I wasn't aware of the fact that nanmax is delivered only in toolboxes. Thanks for that ;-)

Edit: I just found that max also seems to work with nan? I'm confused...

[Guillaume]

Timing demo:

im_raw = randi(1000, 2546, 3318);
im_raw(im_raw > 998) = NaN;
fprintf('original image has %d nans\n\n', sum(isnan(im_raw(:))));
orig_im = im_raw;

tic
im_raw(isnan(im_raw)) = -inf; %replace nans with -Inf so it doesn't mess up |max|
%pad matrix to avoid edge issues:
im_raw2 = ones(size(im_raw) + 18) * -Inf; %assuming a 19x19 neighbourhood
im_raw2(10:end-9, 10:end-9) = im_raw;   %assuming a 19x19 neighbourhood 
[r, c] = find(isinf(im_raw)); %find values to replace
for p = 1:numel(r)
    neighbourhood = im_raw2(r(p):r(p)+18, c(p):c(p)+18);
    im_raw(r(p), c(p)) = max(neighbourhood(:));
end
toc

Displays (on my machine)

original image has 17055 nans

Elapsed time is 0.184098 seconds.

[Antonio]

Thanks a lot!... now it has worked...;-))