在MATLAB中找到矩阵的NAN边界

Question

我有一个非常大的（2019x1678 double）DEM（数字高程模型）文件作为MATLAB中的矩阵。它的边缘包含NAN值。为了说明我的代码中的边缘效果，我必须在DEM周围放置一个1个单元缓冲区（与相邻单元格相同的值）。在存在NAN的地方，我需要找到NAN值的边缘，以构建该缓冲区。我尝试采取两种方法：

在首先我获得了行和列协调所有非nan dem值，并找到每列的第一个和最后一行数字以获取北然后，向南边界，然后找到每行的第一个也是最后一列编号，以获取东部和西边界。我在 sub2ind（）中使用它们来创建我的缓冲区。

[r, c] = find(~isnan(Zb_ext)); %Zb is my DEM matrix
idx = accumarray(c, r, [], @(x) {[min(x) max(x)]});
idx = vertcat(idx{:});

NorthBoundary_row = transpose(idx(:,1)); % the values to fill my buffer with
NorthBoundary_row_ext = transpose(idx(:,1) - 1); % My buffer cells
columnmax = length(NorthBoundary_row);
column1 = min(c);
Boundary_Colu = linspace(column1,column1+columnmax-1,columnmax);
SouthBoundary_row = (transpose(idx(:,2))); % Repeat for south Boundary
SouthBoundary_row_ext = transpose(idx(:,2) + 1);

SouthB_Ind = sub2ind(size(Zb_ext),SouthBoundary_row,Boundary_Colu);
SouthB_Ind_ext = sub2ind(size(Zb_ext),SouthBoundary_row_ext, Boundary_Colu);
NorthB_Ind = sub2ind(size(Zb_ext),NorthBoundary_row, Boundary_Colu);
NorthB_Ind_ext = sub2ind(size(Zb_ext),NorthBoundary_row_ext, Boundary_Colu);

Zb_ext(NorthB_Ind_ext) = Zb_ext(NorthB_Ind);
Zb_ext(SouthB_Ind_ext) = Zb_ext(SouthB_Ind);

% Repeat above for East and West Boundary by reversing the roles of row and 
% column

[r, c] = find(~isnan(Zb_ext));
idx = accumarray(r, c, [], @(x) {[min(x) max(x)]});
idx = vertcat(idx{:});

EastBoundary_colu = transpose(idx(:,1)); % Repeat for east Boundary
EastBoundary_colu_ext = transpose(idx(:,1) - 1); 
row1 = min(r);
rowmax = length(EastBoundary_colu);
Boundary_row = linspace(row1,row1+rowmax-1,rowmax);
WestBoundary_colu = transpose(idx(:,2)); % Repeat for west Boundary
WestBoundary_colu_ext = transpose(idx(:,2) + 1);

EastB_Ind = sub2ind(size(Zb_ext),Boundary_row, EastBoundary_colu);
EastB_Ind_ext = sub2ind(size(Zb_ext),Boundary_row, EastBoundary_colu_ext);
WestB_Ind = sub2ind(size(Zb_ext),Boundary_row, WestBoundary_colu);
WestB_Ind_ext = sub2ind(size(Zb_ext),Boundary_row, WestBoundary_colu_ext);

Zb_ext(NorthB_Ind_ext) = Zb_ext(NorthB_Ind);
Zb_ext(SouthB_Ind_ext) = Zb_ext(SouthB_Ind);
Zb_ext(EastB_Ind_ext) = Zb_ext(EastB_Ind);
Zb_ext(WestB_Ind_ext) = Zb_ext(WestB_Ind);

这在我的小型开发矩阵上效果很好，但在我的全尺寸dem上失败了。我不了解我的代码的行为，但是查看我边界上有空白的数据。我想知道我是否需要更好地控制最大/min行/列值的顺序，尽管在较小的数据集上的测试中，所有这些似乎都按顺序进行。...

我从a 模仿问题扩张方法。但是，当我过渡到完整的数据集时，计算 zbdilated 需要数小时。尽管我的第一种方法不起作用，但它至少在几秒钟内计算。

[m, n] = size(Zb); % 
Zb_ext = nan(size(Zb)+2);
Zb_ext(2:end-1, 2:end-1) = Zb; % pad Zb with zeroes on each side
ZbNANs = ~isnan(Zb_ext);
ZbDilated = zeros(m + 2, n + 2); % this will hold the dilated shape.

for i = 1:(m+2)
    if i == 1 %handling boundary situations during dilation
        i_f = i;
        i_l = i+1;
    elseif i == m+2
        i_f = i-1;
        i_l = i;
    else
        i_f = i-1;
        i_l = i+1;
    end
    for j = 1:(n+2)
        mask = zeros(size(ZbNANs));
        if j == 1 %handling boundary situations again
            j_f = j;
            j_l = j+1;
        elseif j == n+2
            j_f = j-1;
            j_l = j;
        else
            j_f = j-1;
            j_l = j+1;
        end
        
        mask(i_f:i_l, j_f:j_l) = 1; % this places a 3x3 square of 1's around (i, j)
        ZbDilated(i, j) = max(ZbNANs(logical(mask)));
    end
end

Zb_ext(logical(ZbDilated)) = fillmissing(Zb_ext(logical(ZbDilated)),'nearest');

有人对这两个可用的想法有任何想法吗？

这是我首先：

   NaN   NaN     2     5    39    55    44     8   NaN   NaN
   NaN   NaN   NaN     7    33    48    31    66    17   NaN
   NaN   NaN   NaN    28   NaN    89   NaN   NaN   NaN   NaN

以下是与NAN的限制缓冲的矩阵：

   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   NaN     2     5    39    55    44     8   NaN   NaN   NaN
   NaN   NaN   NaN   NaN     7    33    48    31    66    17   NaN   NaN
   NaN   NaN   NaN   NaN    28   NaN    89   NaN   NaN   NaN   NaN   NaN
   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN

这是我在使用 fillmissing 之后要获得的（尽管我注意到填充缓冲区的填充方式有些不规则性。 。

   NaN   NaN     2     2     5    39    55    44     8    17   NaN   NaN
   NaN   NaN     2     2     5    39    55    44     8    17    17   NaN
   NaN   NaN     2     2     7    33    48    31    66    17    17   NaN
   NaN   NaN   NaN     2    28    33    89    31    66    17    17   NaN
   NaN   NaN   NaN     5    28    55    89     8   NaN   NaN   NaN   NaN

​

       0   0   1     1    1    1    1    1   1   1   0   0
       0   0   1     1    1    1    1    1   1   1   1   0
       0   0   1     1    1    1    1    1   1   1   1   0
       0   0   0     1    1    1    1    1   1   1   1   0
       0   0   0     1    1    1    1    1   0   0   0   0

Answer 1

应用3x3扩张的速度将如下。这确实涉及一些大型的中间矩阵，这使其效率较低，例如应用 Imdilate 。

[m, n] = size(Zb); % 
Zb_ext = nan(size(Zb)+2);
Zb_ext(2:end-1, 2:end-1) = Zb; % pad A with zeroes on each side
ZbNANs = ~isnan(Zb_ext);
ZbDilated = ZbNANs; % this will hold the dilated shape.

% up and down neighbors
ZbDilated(2:end, :) = max(ZbDilated(2:end, :), ZbNANs(1:end-1, :));
ZbDilated(1:end-1, :) = max(ZbDilated(1:end-1, :), ZbNANs(2:end, :));

% left and right neighbors
ZbDilated(:, 2:end) = max(ZbDilated(:, 2:end), ZbNANs(:, 1:end-1));
ZbDilated(:, 1:end-1) = max(ZbDilated(:, 1:end-1), ZbNANs(:, 2:end));

% and 4 diagonal neighbors
ZbDilated(2:end, 2:end) = max(ZbDilated(2:end, 2:end), ZbNANs(1:end-1, 1:end-1));
ZbDilated(1:end-1, 2:end) = max(ZbDilated(1:end-1, 2:end), ZbNANs(2:end, 1:end-1));
ZbDilated(2:end, 1:end-1) = max(ZbDilated(2:end, 1:end-1), ZbNANs(1:end-1, 2:end));
ZbDilated(1:end-1, 1:end-1) = max(ZbDilated(1:end-1, 1:end-1), ZbNANs(2:end, 2:end));

这是一种乏味的写作方式，我敢肯定有一个可以写得更短的循环，但是我认为这使意图更加清晰。

[编辑：因为我们在这里处理一个逻辑数组，而不是 max（a，b）我们也可以做 a | B 。 时间是否有任何差异。

---

我不确定 /72710739＃comment128432629_72708722“>@beaker在评论中说是不使用的

mask = zeros(size(ZbNANs));
mask(i_f:i_l, j_f:j_l) = 1; % this places a 3x3 square of 1's around (i, j)
ZbDilated(i, j) = max(ZbNANs(logical(mask)));

，而是要使用

ZbDilated(i, j) = max(ZbNANs(i_f:i_l, j_f:j_l), [], 'all');

[编辑：因为我们在这里处理一个逻辑数组，而不是 max（a，a，[]所有'）我们也可以做任何（a，'all'），应该更快。参见。]