在Matlab中过滤包含NaN的图像？

Posted 2023-03-10

【中文标题】在Matlab中过滤包含NaN的图像？

【英文标题】：Filter image that contains NaNs in Matlab?

【发布时间】：2015-07-02 04:57:35

【问题描述】：

我有一个代表一些数据的二维数组 (doubles)，其中有一堆 NaNs。数据的等高线图如下所示：

所有的空白都是NaNs，灰色菱形是供参考的，填充的轮廓显示了我的数据的形状。当我使用imfilt 过滤数据时，NaNs 会显着咀嚼数据，所以我们最终会得到这样的结果：

您可以看到支持集显着收缩。我不能使用它，因为它已经咀嚼了边缘上一些更有趣的变化（由于我的实验所特有的原因，这些边缘很重要）。

是否有在NaNs 岛内进行过滤的功能，它处理类似于矩形过滤窗口边缘的边缘，而不是仅仅杀死边缘？有点像nanmean 函数，除了卷积图像？

这是我的过滤代码：

filtWidth = 7;
imageFilter=fspecial('gaussian',filtWidth,filtSigma);
%convolve them
dataFiltered = imfilter(rfVals,imageFilter,'symmetric','conv');

以及绘制等高线图的代码：

figure
contourf(dataFiltered); hold on
plot([-850 0 850 0 -850], [0 850 0 -850 0], 'Color', [.7 .7 .7],'LineWidth', 1); %the square (limits are data-specific)
axis equal

在 Mathworks 文件交换 (ndanfilter.m) 中有一些代码接近我想要的，但我相信它只会插入散布在图像内部上的 NaNs ，而不是显示这种岛型效应的数据。

注意：我刚刚找到nanconv.m，它完全符合我的要求，使用非常直观（卷积图像，忽略NaN，很像nanmean 工作）。我已将这部分作为我接受的答案的一部分，并与其他答案的表现进行了比较。

相关问题

Gaussian filtering a image with Nan in Python

【问题讨论】：

数据是什么样的？它只是一个 2D 替身吗？

是的，二维双。更新了问题以明确这一点。

imFilt 最后一行...你的意思是imFilter？

您能否也发布一下您是如何绘制图像的？

我想我有一个解决方案，只需要绘制它XD

【参考方案1】：

我最终使用的技术是 Matlab 文件交换中的函数nanconv.m。它完全符合我的要求：它以忽略 NaNs 的方式运行过滤器，就像 Matlab 的内置函数 nanmean 所做的那样。这很难从函数的文档中破译，这有点神秘。

我是这样使用它的：

filtWidth = 7;
filtSigma = 5;
imageFilter=fspecial('gaussian',filtWidth,filtSigma);
dataFiltered = nanconv(data,imageFilter, 'nanout');

我在下面粘贴nanconv 函数（它被the BSD license 覆盖）。当我有机会时，我会发布图片等，只是想发布我最终为任何对我所做的事情感到好奇的人所做的事情。

与其他答案的比较

使用gnovice's solution，结果在直观上看起来非常好，但边缘上存在一些令人担忧的定量信号。在实践中，超出边缘的图像外推导致我的数据边缘出现许多虚假的高值。

使用krisdestruction's suggestion 用原始数据替换丢失的位，看起来也很不错（尤其是对于非常小的过滤器），但是（根据设计）你最终会在边缘得到未经过滤的数据，这对我来说是个问题应用。

nanconv

function c = nanconv(a, k, varargin)
% NANCONV Convolution in 1D or 2D ignoring NaNs.
%   C = NANCONV(A, K) convolves A and K, correcting for any NaN values
%   in the input vector A. The result is the same size as A (as though you
%   called 'conv' or 'conv2' with the 'same' shape).
%
%   C = NANCONV(A, K, 'param1', 'param2', ...) specifies one or more of the following:
%     'edge'     - Apply edge correction to the output.
%     'noedge'   - Do not apply edge correction to the output (default).
%     'nanout'   - The result C should have NaNs in the same places as A.
%     'nonanout' - The result C should have ignored NaNs removed (default).
%                  Even with this option, C will have NaN values where the
%                  number of consecutive NaNs is too large to ignore.
%     '2d'       - Treat the input vectors as 2D matrices (default).
%     '1d'       - Treat the input vectors as 1D vectors.
%                  This option only matters if 'a' or 'k' is a row vector,
%                  and the other is a column vector. Otherwise, this
%                  option has no effect.
%
%   NANCONV works by running 'conv2' either two or three times. The first
%   time is run on the original input signals A and K, except all the
%   NaN values in A are replaced with zeros. The 'same' input argument is
%   used so the output is the same size as A. The second convolution is
%   done between a matrix the same size as A, except with zeros wherever
%   there is a NaN value in A, and ones everywhere else. The output from
%   the first convolution is normalized by the output from the second 
%   convolution. This corrects for missing (NaN) values in A, but it has
%   the side effect of correcting for edge effects due to the assumption of
%   zero padding during convolution. When the optional 'noedge' parameter
%   is included, the convolution is run a third time, this time on a matrix
%   of all ones the same size as A. The output from this third convolution
%   is used to restore the edge effects. The 'noedge' parameter is enabled
%   by default so that the output from 'nanconv' is identical to the output
%   from 'conv2' when the input argument A has no NaN values.
%
% See also conv, conv2
%
% AUTHOR: Benjamin Kraus ([e-mail], [e-mail])
% Copyright (c) 2013, Benjamin Kraus
% $Id: nanconv.m 4861 2013-05-27 03:16:22Z bkraus $

% Process input arguments
for arg = 1:nargin-2
    switch lower(vararginarg)
        case 'edge'; edge = true; % Apply edge correction
        case 'noedge'; edge = false; % Do not apply edge correction
        case 'same','full','valid'; shape = vararginarg; % Specify shape
        case 'nanout'; nanout = true; % Include original NaNs in the output.
        case 'nonanout'; nanout = false; % Do not include NaNs in the output.
        case '2d','is2d'; is1D = false; % Treat the input as 2D
        case '1d','is1d'; is1D = true; % Treat the input as 1D
    end
end

% Apply default options when necessary.
if(exist('edge','var')~=1); edge = false; end
if(exist('nanout','var')~=1); nanout = false; end
if(exist('is1D','var')~=1); is1D = false; end
if(exist('shape','var')~=1); shape = 'same';
elseif(~strcmp(shape,'same'))
    error([mfilename ':NotImplemented'],'Shape ''%s'' not implemented',shape);
end

% Get the size of 'a' for use later.
sza = size(a);

% If 1D, then convert them both to columns.
% This modification only matters if 'a' or 'k' is a row vector, and the
% other is a column vector. Otherwise, this argument has no effect.
if(is1D);
    if(~isvector(a) || ~isvector(k))
        error('MATLAB:conv:AorBNotVector','A and B must be vectors.');
    end
    a = a(:); k = k(:);
end

% Flat function for comparison.
o = ones(size(a));

% Flat function with NaNs for comparison.
on = ones(size(a));

% Find all the NaNs in the input.
n = isnan(a);

% Replace NaNs with zero, both in 'a' and 'on'.
a(n) = 0;
on(n) = 0;

% Check that the filter does not have NaNs.
if(any(isnan(k)));
    error([mfilename ':NaNinFilter'],'Filter (k) contains NaN values.');
end

% Calculate what a 'flat' function looks like after convolution.
if(any(n(:)) || edge)
    flat = conv2(on,k,shape);
else flat = o;
end

% The line above will automatically include a correction for edge effects,
% so remove that correction if the user does not want it.
if(any(n(:)) && ~edge); flat = flat./conv2(o,k,shape); end

% Do the actual convolution
c = conv2(a,k,shape)./flat;

% If requested, replace output values with NaNs corresponding to input.
if(nanout); c(n) = NaN; end

% If 1D, convert back to the original shape.
if(is1D && sza(1) == 1); c = c.'; end

end

【讨论】：

注意，最近提交给 File Exchange 的实现稍微灵活一些：mathworks.com/matlabcentral/fileexchange/20417-ndnanfilter-m

【参考方案2】：

一种方法是在执行过滤之前使用 scatteredInterpolant（或旧 MATLAB 版本中的 TriScatteredInterp）将 NaN 值替换为最近邻插值，然后再将这些点替换为 NaN 值。这类似于使用'replicate' 参数过滤完整的二维数组，而不是使用'symmetric' 参数作为imfilter 的边界选项（即，您正在复制而不是在锯齿状NaN 边界反射值)。

代码如下所示：

% Make your filter:
filtWidth = 7;
imageFilter = fspecial('gaussian', filtWidth, filtWidth);

% Interpolate new values for Nans:
nanMask = isnan(rfVals);
[r, c] = find(~nanMask);
[rNan, cNan] = find(nanMask);
F = scatteredInterpolant(c, r, rfVals(~nanMask), 'nearest');
interpVals = F(cNan, rNan);
data = rfVals;
data(nanMask) = interpVals;

% Filter the data, replacing Nans afterward:
dataFiltered = imfilter(data, imageFilter, 'replicate', 'conv');
dataFiltered(nanMask) = nan;

【讨论】：

你对使用imfilter(....'replicate')和函数nanconv（mathworks.com/matlabcentral/fileexchange/41961-nanconv）的区别有了解吗？他们都给出了相似的结果，但坦率地说，我并不真正了解其中一个是如何工作的。

看起来很有希望（我的版本我只有TriScatteredInterp，但它的工作原理完全相同）。目前正在将其与nanconv 的性能进行比较，实际上......

@neuronet 我不熟悉nanconv 或者底层代码是使用，所以我还不确定可能有什么区别。

@neuronet：使用您问题中的数据添加一些数字会很好。对于旧版本，我将在提及 TriScatteredInterp 时进行编辑。

@krisdestruction 是的，它是一样的。它只在没有数据点的地方进行插值。今天我一直在玩这个，虽然它看起来很直观，但从数量上看，边缘有点奇怪，因为过滤外插图像会放大边缘的微小变化。我发现的nanconv 函数似乎工作得更好一些，更像是nanmean，因为它只是忽略 nans，很像nanmean 内置函数。这确实是一个非常聪明的功能。很可能会将此作为单独的答案发布，因此忙于工作（我正在分析这样的数据:)）。

【参考方案3】：

好的，不使用你的绘图功能，我仍然可以给你一个解决方案。您要做的是找到所有新的NaN 并将其替换为原始未过滤数据（假设它是正确的）。虽然它没有被过滤，但它比减少轮廓图像的域更好。

% Toy Example Data
rfVals= rand(100,100);
rfVals(1:2,:) = nan;
rfVals(:,1:2) = nan;

% Create and Apply Filter
filtWidth = 3;
imageFilter=fspecial('gaussian',filtWidth,filtWidth);
dataFiltered = imfilter(rfVals,imageFilter,'symmetric','conv');
sum(sum(isnan( dataFiltered ) ) )

% Replace New NaN with Unfiltered Data
newnan = ~isnan( rfVals) & isnan( dataFiltered );
dataFiltered( newnan ) = rfVals( newnan );
sum(sum(isnan( rfVals) ) )
sum(sum(isnan( dataFiltered ) ) )

使用以下代码检测新的NaN。您也可以使用xor 函数。

newnan = ~isnan( rfVals) & isnan( dataFiltered );

然后这一行将dataFiltered 中的索引设置为rfVals 中的值

dataFiltered( newnan ) = rfVals( newnan );

结果

从控制台打印的行和我的代码中，您可以看到 dataFiltered 中 NaN 的数量从 688 减少到 396，rfVals 中 NaN 的数量也是如此。

ans =
   688
ans =
   396
ans =
   396

替代解决方案 1

您还可以通过指定较小的内核并在之后合并它来在边缘附近使用较小的过滤器，但如果您只想要使用最少代码的有效数据，我的主要解决方案将起作用。

替代解决方案 2

另一种方法是用零或您想要的某个常数填充/替换NaN 值，以便它可以工作，然后截断它。但是对于信号处理/过滤，您可能需要我的主要解决方案。

【讨论】：

这会起作用，而且非常好，但会在边缘引入不连续性和异常高的频率。理想情况下，我们可以构建类似于nanmean 的东西，但是一个过滤器在卷积时会进行适当的过滤，但仅限于不是NaN 的数据点。将玩弄你的想法，发布它的样子......

我不想添加零（这就是我使用蒙版所做的），但它会将边缘的值拉低....

@neuronet 是的，这只是一个替代解决方案，我知道你不希望这样。我的主要解决方案会提取您的主要数据中缺少的数据。

@neuronet 另一种替代方法是在边缘附近使用较小尺寸的过滤器。虽然它更准确，但涉及的代码更多。

是的，如果它有帮助并且有效，我当然会接受，无需不断提醒我:) 我正在研究类似nanconvolve 的东西来比较它们。替代解决方案非常好，但奥卡姆的剃刀开始抬头。就像我说的那样，研究看起来最简洁的替代方案，它会卷积但忽略 NaNs。

【参考方案4】：

nanfilter 在过滤时与 nanconv 完全相同，只要过滤器相同。如果您在使用 nanfilter 之前获得了 nan 值，然后将其添加到后过滤矩阵中，只要您使用相同的过滤器，您将获得与使用选项 'nanout' 从 nanconv 获得的结果相同的结果。

【讨论】：

什么是nanfilter？这个页面还没有提到它，我找不到它作为一个内置的matlab函数，谷歌什么也没透露。也许您打算将此添加为评论？