How to pad an irregularly shaped matrix in matlab-CodePudding

I have a matrix with values in the center and NaNs on the border (imagine a matrix representing a watershed which is never square). I need to pad it with one cell to do some component stress calculations. I am trying to avoid using outside libraries from the core Matlab functionality however what i am trying to do is similar to padarray symmetric but for an irregular border:

padarray(Zb,[1 1],'symmetric','both');

For example:

   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

Becomes:

   NaN   NaN     2     2     5    39    55    44     8    8   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    28    28    33    89    31    66    17    17   NaN
   NaN   NaN   NaN    28    28    28    89    89   NaN   NaN   NaN   NaN

(Not sure how to handle convex corners with two adjacent values since I need to control edge effects).

This post follows on an earlier question today in which I was able to extract the locations of these padded cells (buffers) into a dilated logical. However using fillmissing with nearest did not create the effect I expected (what padarray does).

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

I might be able to reverse what I did to find the padcells to find the adjacent values and use those to replace the pad cell NaNs. But I thought I would first see if there was a simpler solution?

CodePudding user response：

You can use two 2D convolutions to achieve this, where conv2 is within the core MATLAB library so nothing external is needed, and it should be fast.

However, you noted this:

Not sure how to handle convex corners with two adjacent values since I need to control edge effects

I've taken the liberty of defining a "sensible" output for convex corners which is to take the average value, because from your example it seems undefined how these cases, and more complicated ones like cell (5,6), should be handled.

I've added detailed comments to the below code for explanation

% Example matrix 
A = [
    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
    ];

% Track the "inner" indices, where values are defined
inner = (~isnan(A));
B = A;                 % Copy A so we don't change it
B(~inner) = 0;         % Replace NaN with 0 so that convolutions work OK

% First dilate the inner region by one element, taking the average of
% neighbours which are up/down/left/right (no diagonals). This is required
% to avoid including interior points (which only touch diagonally) in the
% averaging. These can be considered the "cardinal neighbours"
kernel = [0 1 0 ; 1 0 1; 0 1 0]; % Cardinal directions in 3x3 stencil
s = conv2(B,kernel,'same');      % 2D convolution to get sum of neighbours
n = conv2(inner,kernel,'same');  % 2D convolution to get count of neighbours
s(inner) = 0;                    % Zero out the inner region
s = s./n;                        % Get the mean of neighbours

% Second, dilate the inner region but including the mean from all
% directions. This lets us handle convex corners in the image
s2 = conv2(B,ones(3),'same');     % Sum of neighbours (and self, doesn't matter)
n = conv2(inner,ones(3),'same');  % Count of neighbours (self=0 for dilated elems)
s2 = s2./n;                       % Get the mean of neighbours

% Finally piece together the 3 matrices:
out = s2;                       % Start with outmost dilation inc. corners
out(~isnan(s)) = s(~isnan(s));  % Override with inner dilation for cardinal neighbours
out(inner) = A(inner);          % Override with original inner data

So for this example, the output would be the same as your example output, except for corners as mentioned:

NaN  NaN    2    2     5     39   55   44    8    8   NaN  NaN 
NaN  NaN    2    2     5     39   55   44    8  12.5   17  NaN 
NaN  NaN    2    4.5   7     33   48   31   66   17    17  NaN 
NaN  NaN  NaN   28    28     50   89   60   66   17    17  NaN 
NaN  NaN  NaN   28    28   58.5   89   89  NaN  NaN   NaN  NaN