Help improving speed of this code/loop

Dear all,

I wrote the code below, which does exactly what I need. Basically, I have several (10.000+) 134x134 matrices (images). For each pixel of each matrix I first need to calculate the average and the standard deviation of the surrounding 8 pixels (removing the central pixel from the calculation). Then I create my Lower and Upper boundaries as the Mean-Std and Mean+Std, respectively. Then, if any, I remove the values of the 8 pixels falling outside these boundaries and I recompute the average. Finally, that average value is placed in a new matrix, at the same position of the central pixel.

Although the code does the job, it takes more than 20 seconds to go through the loop. I need to run it more than 10 thousand times, and it takes ages. Would you be able to suggest an alternative, improved code to make it faster? Perhaps removing the loop, or maybe exploiting some built-in or external function?

Thanks a lot for your help!

A=magic(134); % Assuming this is my 134x134 matrix
background=NaN(size(A));
Mtx_Siz=134*134;

xmin=1;
xmax=3;
ymin=1;
ymax=3;

for i=1:Mtx_Siz
    if ymax<135;
        Mask=zeros(size(A));
        Mask(ymin:ymax,xmin:xmax)=1;
        Mask(ymax-1,xmax-1)=0;
        jmask=find(Mask==1);
        B=A(jmask);
        B_Mean=nanmean(B);
        B_std=nanstd(B);
        B_Upper=B_Mean+B_std;
        B_Lower=B_Mean-B_std;
        jout=find(B<B_Lower|B>B_Upper);
        B(jout)=NaN;
        B_Mean2=nanmean(B);
        B_std2=nanstd(B);
        background(ymax-1,xmax-1)=B_Mean2;
        xmin=xmin+1;
        xmax=xmax+1;

        if xmax==135;
            xmin=1;
            xmax=3;
            ymin=ymin+1;
            ymax=ymax+1;

        else
        end
    else
    end
    keep A background Mtx_Siz ymin ymax xmin xmax
end

[Stephen23]

Start by profiling your code: https://www.mathworks.com/help/matlab/matlab_prog/profiling-for-improving-performance.html

Get rid of KEEP (which calls slow CLEAR via slow EVALIN, see the warnings here: https://www.mathworks.com/help/matlab/matlab_prog/techniques-for-improving-performance.html#buwj1nn ). Clearing variables in a loop usually just gets gets in the way of MATLAB taking care of memory efficiently by itself: it is better to write your code so that clearing is not required.

Most/all of those FINDs can be removed as well: logical indexing is simpler and faster.

Perhaps something like BLOCKPROC or CONV2 might help you?

[Simone A.]

Dear Stephen, thanks a lot for your timely and helpful reply.

I was using conv2 (nanconv) to do a simple average of the neighbouring pixels, and it was ultra-fast.

However, I have no idea how to implement the upper and lower limit bit using conv 2. Would you be able to help?

I was able to obtain the simple average with nanconv by:

A=magic(134); % Assuming this is my 134x134 matrix
kernel = ones(3, 3) / 8;
kernel(2, 2) = 0;
background = nanconv(A, kernel, 'same');

Any chance this piece of code could be adapted to include the tasks presented in the original question?

Thanks once again.

Answer 1

[Mathieu NOE]

hello

maybe this runs a bit faster ... my PC is not fast at all , and this loop needs less than 1 s to complete. I am not sure you can do all this in a non for loop code (there could be for loop inside functions anyway).

N = 134;
A=magic(N); % Assuming this is my 134x134 matrix
background=NaN(size(A));
Mtx_Siz=N*N;
xmin=1;
xmax=3;
ymin=1;
ymax=3;

% some preparation work to ease the for loop job
[r,c] = ind2sub(size(A),1:Mtx_Siz); % create r,c indexes for all pixels
% remove first & last row and col indexes from the list
idr = (r>r(1) & r<r(end));
idc = (c>c(1) & c<c(end)); 
id = idr & idc;
r = r(id);
c = c(id);

tic

for i=1:numel(r)
        ri = r(i)-1:r(i)+1;
        ci = c(i)-1:c(i)+1;
        B=A(ri,ci);
        B = B(:); % vector of lenght = 9 , central pixel is at 5th position
        B(5) = [];% remove central pixel value
        B_Mean=mean(B,'omitnan');
        B_std=std(B,'omitnan');
        B_Upper=B_Mean+B_std;
        B_Lower=B_Mean-B_std;
        jout=(B<B_Lower|B>B_Upper);
        B(jout)=NaN;
        B_Mean2=mean(B(:),'omitnan');
        B_std2=std(B(:),'omitnan');
        background(ri,ci)=B_Mean2;

end

toc

Elapsed time is 0.404206 seconds.

[Mathieu NOE]

hello gain

did my suggestion bring some benefit on your side ?