Help increasing the speed of a piece of code

2 views (last 30 days)
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

Accepted Answer

Rik
Rik on 15 Apr 2023
Edited: Rik on 15 Apr 2023
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.
  3 Comments
Rik
Rik on 15 Apr 2023
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.
Simone A. on 16 Apr 2023
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.

Sign in to comment.

More Answers (1)

per isakson
per isakson on 15 Apr 2023
Edited: per isakson on 15 Apr 2023
"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.
  1 Comment
Simone A.
Simone A. on 16 Apr 2023
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

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!