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

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);

Answer 1

[William Rose]

@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.

[DGM]

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]

@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]

@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.]

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]

@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.

Answer 2

[Sulaymon Eshkabilov]

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

Answer 3

[DGM]

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?