Help increasing the speed of a piece of code

Dear all,

Plain summary: I have a temperature matrix and I need to find the average temperature of the 8 pixels surrounding each pixel.

I need to subdivide my 134x134 matrices into 3x3 sub matrices. For each submatrix I need to compute the mean and the standard deviation, excluding the central pixel from the calculation. Then I need to find the values exceding the upper and lower boundaries, namely the mean + standard deviation and the mean - standard deviation, respectively. Then, I need to recalculate the mean excluding the values falling outside the lower and upper boundaries and assign the mean value to the central pixel location (in a new matrix).

The output will be a 134x134 matrix with every value (pixel) corresponding to the average of its surrounding pixels with the "outliers" excluded from the average.

I wrote a code (loop) to do the job, but it was taking ~20 seconds (link to old question and previous code: https://it.mathworks.com/matlabcentral/answers/1947148-help-improving-speed-of-this-code-loop?s_tid=srchtitle ). Someone helped me to improve the original code (see new code below), but it still takes ~2 seconds.

As I need to run it 10k+ times it will take ages to complete. I wonder if you could kindly suggest a faster code to achieve the same result?

Thanks a lot!

A=magic(134); % Assuming this is my 134x134 Temperature matrix

T_background = blockproc(A,[1,1],@BlkFun, 'BorderSize',[1,1],...
    'TrimBorder',false, 'PadMethod',NaN);

%where BlkFun=

function v = BlkFun(s)
px = s.data;
if isequal(size(px),[3,3])
    px(5) = [];
    av = nanmean(px);
    sd = nanstd(px);
    ix = px>=(av-sd) & px<=(av+sd);
    v = nanmean(px(ix));
else
    v = NaN;
end
end

Answer 1

[Rik]

It sounds like a convolution would be perfect here.

kernel=ones(3,3);
kernel(2,2)=0;
kernel=kernel/sum(kernel(:));

That will find the average of the surrounding pixels.

Another strategy would be to create a 3D array with all shifts. Then you can use the functions already in your function. Instead of removing elements (or indexing) you should mark invalid values with NaN. Memory limits should not be an issue since it is only 134x134x8.

[Simone A.]

Hi @Rik thanks for your reply. Yes, conv2 (or nanconv) does the job for the avg and it's extremely quick but I need to do the other steps(std, upeer and lower bounds, etc.). I have no I dea how to do that.

About "create a 3D array with all shifts", I have no idea how to do that. Would you be able to adapt my code to it? That would be great!

Thanks a lot for your time.

[Rik]

A=magic(4); % Assuming this is my 134x134 Temperature matrix

% pad the array with NaNs to make circshift easier
A_pad = [...
    NaN(1,2+size(A,2));...
    NaN(size(A,1),1) A NaN(size(A,1),1);...
    NaN(1,2+size(A,2))];
A_3D = repmat(A_pad,[1 1 8]); % Extend to 3D
k = 0;
for row_shift = [-1 0 1]
    for col_shift = [-1 0 1]
        % skip the center pixel
        if row_shift==0 && col_shift==0,continue,end
        k = k+1;
        A_3D(:,:,k) = circshift(A_pad,[row_shift col_shift]);
    end
end
%A_3D([1 end],[1 end],:) = []; % remove NaN padding to get back to the original size
A_3D([1 end],:,:) = [];A_3D(:,[1 end],:) = []; % Matlab complains if you do this in one go
A_3D

A_3D = 
A_3D(:,:,1) =

    11    10     8   NaN
     7     6    12   NaN
    14    15     1   NaN
   NaN   NaN   NaN   NaN


A_3D(:,:,2) =

     5    11    10     8
     9     7     6    12
     4    14    15     1
   NaN   NaN   NaN   NaN


A_3D(:,:,3) =

   NaN     5    11    10
   NaN     9     7     6
   NaN     4    14    15
   NaN   NaN   NaN   NaN


A_3D(:,:,4) =

     2     3    13   NaN
    11    10     8   NaN
     7     6    12   NaN
    14    15     1   NaN


A_3D(:,:,5) =

   NaN    16     2     3
   NaN     5    11    10
   NaN     9     7     6
   NaN     4    14    15


A_3D(:,:,6) =

   NaN   NaN   NaN   NaN
     2     3    13   NaN
    11    10     8   NaN
     7     6    12   NaN


A_3D(:,:,7) =

   NaN   NaN   NaN   NaN
    16     2     3    13
     5    11    10     8
     9     7     6    12


A_3D(:,:,8) =

   NaN   NaN   NaN   NaN
   NaN    16     2     3
   NaN     5    11    10
   NaN     9     7     6

Now you can use this array to calculate the SD and mean to your heart's content:

mean(A_3D,3,'omitnan')

ans = 4×4

6.0000 9.0000 8.8000 7.0000 9.0000 7.2500 7.7500 8.8000 8.2000 9.2500 9.7500 8.0000 10.0000 8.2000 8.0000 11.0000

[Simone A.]

Dear @Rik thanks a lot for solving my issue! This is exactly what I was after. The final code (attached below for future reference) is ultra fast and does the job!

tic
A=magic(134);
A_pad = [...
    NaN(1,2+size(A,2));...
    NaN(size(A,1),1) A NaN(size(A,1),1);...
    NaN(1,2+size(A,2))];
A_3D = repmat(A_pad,[1 1 8]); % Extend to 3D
k = 0;
for row_shift = [-1 0 1]
    for col_shift = [-1 0 1]
        % skip the center pixel
        if row_shift==0 && col_shift==0,continue,end
        k = k+1;
        A_3D(:,:,k) = circshift(A_pad,[row_shift col_shift]);
    end
end
A_3D([1 end],:,:) = [];A_3D(:,[1 end],:) = []; % Matlab complains if you do this in one go

av=mean(A_3D,3,'omitnan');
sd = std(A_3D,[],3,"omitnan");
UpperL=av+sd;
LowerL=av-sd;

A_3D_Temp=A_3D;
mask_temp = (A_3D_Temp >= UpperL | A_3D_Temp <= LowerL);
A_3D_Temp(mask_temp) = NaN;
background    = mean(A_3D_Temp,3,'omitnan');

toc

% Elapsed time is 0.011887 seconds.

Answer 2

[per isakson]

"I need to do the other steps". Removing NaN before calculating the mean cuts the elapse time in half (on my PC with R2018b)

tic
A=magic(134); % Assuming this is my 134x134 Temperature matrix
T_background = blockproc(A,[1,1],@BlkFun, 'BorderSize',[1,1],...
    'TrimBorder',false, 'PadMethod',NaN);
toc

Elapsed time is 0.404461 seconds.

function v = BlkFun(s)
px = s.data;
if isequal(size(px),[3,3])
    px(isnan(px)|[0,0,0;0,1,0;0,0,0]) = [];
    av = mean(px);
    sd = std(px);
    ix = px>=(av-sd) & px<=(av+sd);
    v  = mean(px(ix));
else
    v = NaN;
end
end

and hopefully it returns the same result.

[Simone A.]

Dear @per isakson thanks a lot for taking the time to and for improving my code. That adjustments cut the time in half indeed. However, I found a way faster alternative by adapting the code suggested by @Rik.

Find the final (and very fast) code above (if needed for future reference).

Thanks once again