Moving Standard Deviation issues with NaN values (i.e., Stdfilt). Alternatives and workarounds?
Show older comments
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);
Accepted Answer
William Rose
on 19 Dec 2023
Edited: William Rose
on 19 Dec 2023
[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 Comments
DGM
on 19 Dec 2023
This will still throughput NaNs if there are enough NaNs in a region to span the filter. For the example I gave with a 10x10 NaN block, the output will have a 6x6 NaN block (10 - 2*floor(5/2) = 6).
It's arguable that's appropriate, but the intent seems to be avoiding NaN throughput -- though stdfilt() already does prevent NaN throughput, so I'm not sure of anything.
William Rose
on 20 Dec 2023
@DGM, you are right (as usual): if the image chunk A is all NaNs, std(A(~isnan(A))) will return NaN.
@Simone A., if that could really happen with your images, the you must decide what value you want to report in such cases, and replace the line
img2std(i,j)=std(A(~isnan(A))); % stdev(A without NaNs)
with something such as
if sum(~isnan(A))>0
img2std(i,j)=std(A(~isnan(A))); % stdev(A without NaNs)
else
img2std(i,j)=0; % 0 if A is all NaNs
end
Good luck.
William Rose
on 20 Dec 2023
@Simone A., in case it is not obvious, img2std, which is displayed as the right hand panel in my answer, does not have any NaN values, even though the source image (in the middle, above) has about 845 NaNs randomly scattered. Even if you make the source image 50% NaNs, there are no NaNs in the output image, and the output still looks decent - see below.
The key here is that the NaNs are randomly scattered, so the chance that all 24 pixels are NaN in any 5x5-1 chunk is very small.
Simone A.
on 20 Dec 2023
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
on 20 Dec 2023
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.
More Answers (2)
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
%% (2) Substitute NAN with some value
SOME_Value = 2;
A1 = A;
A1(isnan(A1))= SOME_Value
%% (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(isnan(A2)) = A2(find(isnan(A2))+1) % In case, there are two nans one after another
%% (3) Skip or omit NAN while computing std()
STD_A = std(A, "omitnan")
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
% 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
% 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
% see the output region local to the defect
B(16:35,16:35)
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
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?
Categories
Find more on Region and Image Properties in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
