Moving Standard Deviation issues with NaN values (i.e., Stdfilt). Alternatives and workarounds?

8 次查看（过去 30 天）
Simone A.
Simone A. 2023-12-19
评论： William Rose 2023-12-20
Hi All,
I need to compute a 5x5 moving standard deviation in a matrix I. Also, I need to exclude the central pixel of my moving window from the std computation (see the code below).
The issue is when I have some NaN values, as the stdfilt output returns a sqaure of NaN.
Can you suggest a solution?
Thanks a lot in advance!
(Please see attached my code below)
I = magic(130);
Stdlfilt_neighborhoodSize = [5, 5];
mask_Stdlfilt = true(Stdlfilt_neighborhoodSize);
mask_Stdlfilt(ceil(Stdlfilt_neighborhoodSize(1)/2), ceil(Stdlfilt_neighborhoodSize(2)/2)) = false;
Std_I = stdfilt(I, mask_Stdlfilt);

请先登录，再回答此问题。

采纳的回答

William Rose
William Rose 2023-12-19
编辑：William Rose 2023-12-19
@Simone A.,
[Edit: I attached the script twice. I am removing the earlier version.]
See attached script. Key steps, which occur inside a loop, are shown below. In the code below, img2 is the image, with NaNs. img2std is the array of standard deviations of img2. sn is the number of points in the st.dev mask. shw is the half-width of the st.dev. mask, rounded down. When the mask is 5x5, sn=25 and shw=2.
A=img2(i-shw:i+shw,j-shw:j+shw); % extract the chunk of the image
A=reshape(A,sn,1); % make A a vector
A((sn+1)/2)=[]; % delete the central element
img2std(i,j)=std(A(~isnan(A))); % stdev(A without NaNs)
The code above computes the st.dev. of A, with the central element and any NaNs removed.
This script generates and displays three images, shown below: the original "magic(130)" image, the original plus NaNs at a random 5% of the pixels, and the st.dev. across the image.
  5 个评论
显示 3更早的评论
Simone A.
Simone A. 2023-12-20
Hi @William Rose Thanks a lot for your help and support! Indeed, your code does exactly what I needed. In the unlikely scenario I should encounter that many NaNs all clustered together, i'll interpolate it.
However, I came up with 2 other alternatives, as I also need to maximise the execution speed.
From testing, the faster code is Method 2 (I set yours as Method 3). Please find the 3 codes below for future reference. Yet, I noticed that the faster method returns slightly different results as compared to the other 2. Discrepancies are neglegible, in the order of >10^-12, but I wonder why they are...
That said, find the code attached below (Sorry it's a bit mess but I had to rush):
% Execution time Method 1 : 0.179709 seconds.
% Execution time Method 2 : 0.010063 seconds.
% Execution time Method 3 : 0.177544 seconds.
clc
close
clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Method 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
II = magic(130);
II(67,67:68) = NaN;
neighborhoodSize = [5, 5]; % Define the neighborhood size
mask = ones(neighborhoodSize); % Create a custom mask to exclude the central pixel from the neighborhood
mask(ceil(neighborhoodSize(1)/2), ceil(neighborhoodSize(2)/2)) = 0;
stdImage = NaN(size(II)); % Initialize the output image
tic
for i = 1:size(II, 1) % Loop through each pixel in the image
for j = 1:size(II, 2)
if i > neighborhoodSize(1)/2 && i <= size(II, 1) - neighborhoodSize(1)/2 && ...
j > neighborhoodSize(2)/2 && j <= size(II, 2) - neighborhoodSize(2)/2 % Check if the current pixel is not at the border
localNeighborhood = II(i-2:i+2, j-2:j+2); % Extract the local neighborhood
localNeighborhood(3,3) = NaN;
stdValue = std(localNeighborhood(~isnan(localNeighborhood))); % Apply the custom std function excluding NaN values
stdImage(i, j) = stdValue; % Assign the computed standard deviation to the corresponding pixel
end
end
end
toc
Method_1 = stdImage;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Method 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
A_pad = [...
NaN(1,2+size(II,2));...
NaN(size(II,1),1) II NaN(size(II,1),1);...
NaN(1,2+size(II,2))];
A_3D = repmat(A_pad,[1 1 24]); % Extend to 3D
k = 0;
for row_shift = [-2 -1 0 1 2]
for col_shift = [-2 -1 0 1 2]
if row_shift==0 && col_shift==0,continue,end % skip the center pixel
k = k+1;
A_3D(:,:,k) = circshift(A_pad,[row_shift col_shift]);
end
end
A_3D([1 end],:,:) = [];
A_3D(:,[1 end],:) = [];
Std_Neighb=std(A_3D,[],3,'omitnan');
Std_Neighb(:, 1:2) = NaN; Std_Neighb(:, end-1:end) = NaN;
Std_Neighb(1:2, :) = NaN; Std_Neighb(end-1:end, :) = NaN;
toc
Method_2 = Std_Neighb;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Method 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
stdSz=5; % size of standard deviation area (odd)
shw=(stdSz(1)-1)/2; % half-width, rounded down
sn=stdSz^2; % number of points in stdev area
[n,m]=size(II);
img2std=zeros(n,m); % allocate array for st.dev. map
for i=shw+1:n-shw
for j=shw+1:m-shw
% Ranges for i, j are chosen to prevent out-of-range
% errors at the edges
A=II(i-shw:i+shw,j-shw:j+shw);
A=reshape(A,sn,1); % make A a vector
A((sn+1)/2)=[]; % delete the central element
img2std(i,j)=std(A(~isnan(A))); % stdev(A without NaNs)
end
end
toc
Method_3 = img2std;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(1,3,1)
imagesc(Method_1 - Method_2)
colormap gray
axis square
colorbar
title ('Method 1 - Method 2')
subplot(1,3,2)
imagesc(Method_1 - Method_3)
colormap gray
axis square
colorbar
title ('Method 1 - Method 3')
subplot(1,3,3)
imagesc(Method_2 - Method_3)
colormap gray
axis square
colorbar
title ('Method 2 - Method 3')
impixelinfo
linkaxes
William Rose
William Rose 2023-12-20
@Simone A.,
You are welcome.
I don't know why Method 2 gives st.devs. that are sometimes a tiny bit different than the values reported by the Methods 1 and 3. The max StDev is 7714, the mean StDev is 699, so differences of ~10^-12 are irrelevant.
Method 2 is interesting. The idea is: 1. Make an augmented image by adding a 1-pixel-wide strip of NaNs all around the original image. 2. Make a 3D array which is a 24-layer stack of copies of the augmented image. 3. Do a circular shift of each of the 24 layers, so that a drill-down through the stack will hit the 24 points of the 5x5 mask (minus the central point). 4. Compute the StDev of each 24-layer stack of pixels.

请先登录，再进行评论。

更多回答（2 个）

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2023-12-19
Ran in:
NaN values can be substituted with (1) "some value" or (2) next to it value or (3) skip nan values element, e.g.:
A = randi([-3, 3], 10, 5);
B = 0./A;
A(isnan(B))= NaN
A = 10×5
-1 -1 -2 2 -3 3 3 NaN -3 -3 1 -2 1 3 -2 NaN 2 -2 NaN 1 -3 3 -1 -1 -3 -1 -1 -3 NaN NaN 2 3 2 2 -2 -2 -2 1 2 -3 -3 3 NaN -1 2 -2 -2 3 1 -2
%% (2) Substitute NAN with some value
SOME_Value = 2;
A1 = A;
A1(isnan(A1))= SOME_Value
A1 = 10×5
-1 -1 -2 2 -3 3 3 2 -3 -3 1 -2 1 3 -2 2 2 -2 2 1 -3 3 -1 -1 -3 -1 -1 -3 2 2 2 3 2 2 -2 -2 -2 1 2 -3 -3 3 2 -1 2 -2 -2 3 1 -2
%% (2) with the next value
IDX_nan = find(isnan(A));
IDX_NOT_nan = find(~isnan(A));
A2 = A;
A2(IDX_nan) = A2(IDX_nan+1)
A2 = 10×5
-1 -1 -2 2 -3 3 3 1 -3 -3 1 -2 1 3 -2 -3 2 -2 -1 1 -3 3 -1 -1 -3 -1 -1 -3 2 -2 2 3 2 2 -2 -2 -2 1 2 -3 -3 3 3 -1 2 -2 -2 3 1 -2
A2(isnan(A2)) = A2(find(isnan(A2))+1) % In case, there are two nans one after another
A2 = 10×5
-1 -1 -2 2 -3 3 3 1 -3 -3 1 -2 1 3 -2 -3 2 -2 -1 1 -3 3 -1 -1 -3 -1 -1 -3 2 -2 2 3 2 2 -2 -2 -2 1 2 -3 -3 3 3 -1 2 -2 -2 3 1 -2
%% (3) Skip or omit NAN while computing std()
STD_A = std(A, "omitnan")
STD_A = 1×5
2.1794 2.3664 2.1671 2.0659 1.8708
DGM
DGM 2023-12-19
编辑：DGM 2023-12-19
Ran in:
IPT stdfilt() doesn't behave as you describe in any version I've tested.
% a float array with a 10x10 NaN block
A = magic(130);
A(21:30,21:30) = NaN;
nansin = nnz(isnan(A)) % 10x10
nansin = 100
% a 5x5 box filter with a 1px hole
fsz = [5, 5];
fk = true(fsz);
fk(ceil(fsz(1)/2), ceil(fsz(2)/2)) = false;
% the filtered result has no NaNs
B = stdfilt(A, fk);
nansout = nnz(isnan(B)) % no nans
nansout = 0
% the 10x10 NaN block turns into a 14x14 block of zeros
% due to the original block size and the filter radius
% i.e. 10 + 2*floor(5/2) = 14
sqrt(nnz(B == 0)) % block width
ans = 14
% see the output region local to the defect
B(16:35,16:35)
ans = 20×20
1.0e+03 * 1.7666 2.0144 2.0144 1.7666 1.3462 0.8156 0.1387 0.1301 0.1200 0.1147 0.1161 0.1238 0.1338 0.1387 0.1387 5.1894 6.2974 6.2974 5.1894 0.1387 2.0144 2.0144 1.7666 1.3462 0.8156 0.1387 0.1338 0.1238 0.1161 0.1147 0.1200 0.1301 0.1387 0.1387 0.1387 5.1894 6.2974 6.2974 5.1894 0.1387 2.0144 1.7666 1.3462 0.8156 0.1387 0.1387 0.1301 0.1200 0.1147 0.1161 0.1238 0.1338 0.1387 0.1387 0.1387 5.1894 6.2974 6.2974 5.1894 0.1387 1.7666 1.3462 0.8156 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6.2974 5.1894 0.1387 1.3462 0.8156 0.1387 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6.2974 5.1894 0.1387 0.8156 0.1387 0.1387 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6.2974 5.1894 0.1387 0.1387 0.1387 0.1387 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6.2974 5.1894 0.1387 0.1387 0.1387 0.1387 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6.2974 5.1894 0.1387 0.1387 0.1387 0.1387 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6.2974 5.1894 0.1387 0.1387 0.1387 0.1338 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6.2974 5.1894 0.1387
Internally, imfilter() is used to do the work, and while the intermediate results from imfilter() will have NaNs, the final output winds up going through
% where do the NaNs turn into zeros?
intermediateresult = [0 100 225 NaN];
output = sqrt(max(intermediateresult,0)) % this is why
output = 1×4
0 10 15 0
which will return 0 for NaNs in the intermediate result.
Either I'm misreading something, or there's something missing from this picture. Is there some other way in which you're using stdfilt() which is producing different results? Are your NaNs sparse, or are they located in large contiguous regions which can span the filter window? Is avoiding NaN throughput the goal, or is it important to instead preserve existing NaNs?

类别

在 Help Center 和 File Exchange 中查找有关 Orange 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by