您可以首先使用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，并且这些值包含在中位数中。另外，如果元素个数中位数是偶数，那么我就简单的选择歧义右边的元素作为最终输出。
这可能不是您特别想要的。更有效的方法是单独排序所有列，同时保持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，您可以在其中具体确定遇到 nans 时要做什么。如果将标志设置为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 解决方案快得多。如果你要多次调用这段代码，可以考虑使用他的实现：Efficient Implementation of `im2col` and `col2im`
计时测试
为了确定建议的方法是否更快，我将使用timeit 执行时序测试。首先，我将创建一个设置公共变量的函数，创建两个函数，其中第一个是使用nlfilter 的原始方法，第二个是使用建议的方法。我将使用 'omitnan' 的方法，因为它会产生你想要的结果。
这是我写的函数。我已经生成了一个 300 x 300 的输入，就像你设置的那样，它包含 0 到 1 之间的所有随机数。我已经做到了，这个输入中大约 20% 的数字有nan。我还设置了您使用nlfilter 使用的匿名函数来过滤没有nans 的中位数以及邻域大小，即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
如您所见，新方法的运行速度大约比 nlfilter 快 (1.033838 / 0.038697) = 26.7162 倍。还不错！