作为我项目的一部分，我需要使用在rxr窗口上执行中值过滤的代码，并忽略nan值 .

我目前使用MATLAB的nlfilter函数 . 问题是它非常慢：300x300示例需要将近5秒，而MATLAB的medfilt2需要0.2秒 . 有没有人有更高效和优雅的解决方案？

注意：在我的情况下，在边界上的行为并不重要 . 在这个例子中，nlfilter自动用零填充数组，但其他解决方案，如边界重复也是可以的 .

code example:

%Initialize input
r = 3; % window size is 3x3
I = [9,1,6,10,1,5,4;2,4,3,8,8,NaN,5;4,5,8,6,2,NaN,3;5,NaN,6,4,NaN,4,9;3,1,10,9,4,3,2;10,9,10,10,6,NaN,5;10,9,4,1,2,7,2];

%perform median filter on rxr window, igonre nans
f = @(A)median(A(~isnan(A)));
filteredRes = nlfilter(I, [r r], f);
filteredRes(nanMask) = nan;

Expected Results

在过滤之前：

I =
 9     1     6    10     1     5     4
 2     4     3     8     8   NaN     5
 4     5     8     6     2   NaN     3
 5   NaN     6     4   NaN     4     9
 3     1    10     9     4     3     2
10     9    10    10     6   NaN     5
10     9     4     1     2     7     2

过滤后：

filteredRes =
     0    2.0000    3.0000    3.0000    3.0000    2.5000         0
2.0000    4.0000    6.0000    6.0000    6.0000       NaN    3.0000
3.0000    4.5000    5.5000    6.0000    5.0000       NaN    3.0000
2.0000       NaN    6.0000    6.0000       NaN    3.0000    2.5000
2.0000    7.5000    9.0000    7.5000    4.0000    4.0000    2.5000
3.0000    9.0000    9.0000    6.0000    5.0000       NaN    2.0000
     0    9.0000    4.0000    2.0000    1.5000    2.0000         0

谢谢！

1 回答

您可以首先使用padarray填充图像，在每侧填充 floor(r/2) 像素，然后使用im2col重新构建填充图像，以便像素的每个邻域放置在单独的列中 . 接下来，您需要首先将所有 nan 值设置为虚拟值，这样您就不会干扰中值计算...也许就像零一样 . 之后，找到每列的中位数，然后重新塑造成适当大小的图像 .

这样的事情应该有效：

r = 3;
nanMask = isnan(I); % Define nan mask
Ic = I;
Ic(nanMask) = 0; % Create new copy of image and set nan elements to zero
IP = padarray(Ic, floor([r/2 r/2]), 'both'); % Pad image
IPc = im2col(IP, [r r], 'sliding'); % Transform into columns
out = reshape(median(IPc, 1), size(I,1), size(I,2)); % Find median of each column and reshape back
out(nanMask) = nan; % Set nan elements back

我们得到：

>> out

out =

     0     2     3     3     1     1     0
     2     4     6     6     5   NaN     0
     2     4     5     6     4   NaN     0
     1   NaN     6     6   NaN     3     2
     1     6     9     6     4     4     2
     3     9     9     6     4   NaN     2
     0     9     4     2     1     2     0

使用上述方法，与预期结果略有不同的是，我们将所有 nan 值设置为0，并且这些值包含在中位数中 . 另外，如果元素的数量在中位数中是偶数，那么我只选择模糊度右侧的元素作为最终输出 .

这可能不是你想要的具体 . 更有效的方法是单独 sort 所有列，同时保持 nan 值不变，然后确定每列有效的最后一个元素，并确定每个元素，确定中途点的位置，并从排序中选择那些列 . 使用 sort 的好处是 nan 值被推向数组的末尾 .

像这样的东西可以工作：

r = 3;
nanMask = isnan(I); % Define nan mask
IP = padarray(I, floor([r/2 r/2]), 'both'); % Pad image
IPc = im2col(IP, [r r], 'sliding'); % Transform into columns
IPc = sort(IPc, 1, 'ascend'); % Sort the columns
[~,ind] = max(isnan(IPc), [], 1); % For each column, find the last valid number
ind(ind == 1) = r*r; % Handles the case when there are all valid numbers per column
ind = ceil(ind / 2); % Find the halfway point
out = reshape(IPc(sub2ind(size(IPc), ind, 1:size(IPc,2))), size(I,1), size(I,2)); % Find median of each column and reshape back
out(nanMask) = nan; % Set nan elements back

我们现在得到：

>> out

out =

     0     2     3     3     5     4     0
     2     4     6     6     6   NaN     3
     4     5     6     6     6   NaN     3
     3   NaN     6     6   NaN     3     3
     3     9     9     9     4     4     3
     3     9     9     6     6   NaN     2
     0     9     4     2     2     2     0

次要注意事项

最新版本的MATLAB有一个名为 nanflag 的可选第三个输入，您可以在其中明确确定遇到 nan 时要执行的操作 . 如果将标志设置为 omitnan ，则将忽略其计算中的所有 nan 元素，其中默认值为 includenan ，您不必指定第三个参数 . 如果您在中值过滤器调用中指定 omitnan 以及在第一步中跳过 nan 值的设置为0部分，您将从 nlfilter 的输出中得到您想要的内容：

r = 3;
nanMask = isnan(I); % Define nan mask
IP = padarray(I, floor([r/2 r/2]), 'both'); % Pad image
IPc = im2col(IP, [r r], 'sliding'); % Transform into columns
out = reshape(median(IPc, 1, 'omitnan'), size(I,1), size(I,2)); % Find median of each column and reshape back
out(nanMask) = nan; % Set nan elements back

我们得到：

>> out

out =

         0    2.0000    3.0000    3.0000    3.0000    2.5000         0
    2.0000    4.0000    6.0000    6.0000    6.0000       NaN    3.0000
    3.0000    4.5000    5.5000    6.0000    5.0000       NaN    3.0000
    2.0000       NaN    6.0000    6.0000       NaN    3.0000    2.5000
    2.0000    7.5000    9.0000    7.5000    4.0000    4.0000    2.5000
    3.0000    9.0000    9.0000    6.0000    5.0000       NaN    2.0000
         0    9.0000    4.0000    2.0000    1.5000    2.0000         0

更高效的im2col解决方案

用户Divakar实现了更快的版本 im2col ，他已经对其进行了基准测试，并且显示出比MATLAB提供的 im2col 解决方案快得多's image processing toolbox. If you'将多次调用此代码，请考虑使用他的实现：Efficient Implementation of im2col and col2im

时间测试

为了确定提议的方法是否更快，我将使用timeit执行时序测试 . 首先，我将创建一个设置公共变量的函数，创建两个函数，其中第一个是 nlfilter 的原始方法，第二个方法是建议的方法 . 我将使用 'omitnan' 使用该方法，因为它可以产生您想要的结果 .

这里's the function I wrote. I' ve生成了一个300 x 300的输入，就像你如何制作它一样，这个输入中约有20％的数字有 nan . 我还设置了你正在使用 nlfilter 的匿名函数来过滤没有 nan 的中位数以及邻域大小，即3 x 3.然后我在这段代码中定义了两个函数 - 代码执行的原始方法使用 nlfilter 进行过滤以及我在上面使用 omitnan 选项提出的建议：

function time_nan

% Initial setup
rng(112234);
I = rand(300,300);
I(I < 0.2) = nan; % Modify approximately 20% of the values in the input with nan
r = 3; % Median filter of size 3
nanMask = isnan(I); % Determine which locations are nan
f = @(A)median(A(~isnan(A))); % Handle to function used by nlfilter

    function original_method
        filteredRes = nlfilter(I, [r r], f);
        filteredRes(nanMask) = nan;
    end

    function new_method
        IP = padarray(I, floor([r/2 r/2]), 'both'); % Pad image
        IPc = im2col(IP, [r r], 'sliding'); % Transform into columns
        out = reshape(median(IPc, 1, 'omitnan'), size(I,1), size(I,2)); % Find median of each column and reshape back
        out(nanMask) = nan; % Set nan elements back
    end

t1 = timeit(@original_method);
t2 = timeit(@new_method);

fprintf('Average time for original method: %f seconds\n', t1);
fprintf('Average time for new method: %f seconds\n', t2);

end

我目前的机器是HP ZBook G5，配备16 GB RAM和Intel Core i7 @ 2.80 GHz . 当您运行此代码时，我得到以下结果：

Average time for original method: 1.033838 seconds
Average time for new method: 0.038697 seconds

如您所见，新方法的运行速度约为（1.033838 / 0.038697）=比 nlfilter 快26.7162倍 . 不错！

回复于 2024-04-19T05:42:56+08:00

2D中值滤波器，忽略nan值

1 回答

次要注意事项

更高效的im2col解决方案

时间测试

相关问题