How to remove extrapolated regions of a grid created with scatteredInterpolant?

I have a cloud of (evenly or not evenly spaced) points defining the bathymetry of a harbour. The coordinates of the points are given in three column arrays in datapoints.mat .

For computational purposes, I need to resample them over a grid with a used-defined space discretization (say, 5 m). scatteredInterpolant seems to do the job quite well for grid points within the boundaries of the original cloud; however, I still need the grid points falling outside the limits of the original dataset to be NaNs. I tried to put the 'ExtrapolationMethod' option of scatteredInterpolant to 'none', but this works only for grid points falling outside the convex hull of the original dataset.

dx = 5;   %grid step

load('datapoints.mat');

F = scatteredInterpolant(x,y,z,'natural','none'); % create the interpolant

[X,Y] = meshgrid(min(x):dx:max(x),min(y):dx:max(y)); % create a grid over which to interpolate

Z = F(X,Y); % calculate the interpolate values over (X,Y)

subplot(1,2,1)
scatter3(x,y,z,4,z,'.'); view(2); grid on;
axis equal tight; xlabel('x (m)'); ylabel('y (m)');
title('original dataset');
subplot(1,2,2)
mesh(X,Y,Z); view(2); grid on;
title('interpolated grid');
axis equal tight; xlabel('x (m)'); yticklabels('');

To solve this issue, I thought of a (cumbersome) solution in which I evaluate the distance of any given interpolated grid point to the closest cloud point; if this distance is greater than a given threshold, that grid point is set to NaN. It can be computationally heavy, however, if the grid has a lot of points. Is there a more optimal way (maybe already implemented in Matlab) that I am not aware of?

Answer 1

[Steve Eddins]

You could try using an alpha shape to postprocess (mask) Z.

Original code

dx = 5;   %grid step
load('datapoints.mat');
F = scatteredInterpolant(x,y,z,'natural','none'); % create the interpolant
[X,Y] = meshgrid(min(x):dx:max(x),min(y):dx:max(y)); % create a grid over which to interpolate
Z = F(X,Y); % calculate the interpolate values over (X,Y)
subplot(1,2,1)
scatter3(x,y,z,4,z,'.'); view(2); grid on;
axis equal tight; xlabel('x (m)'); ylabel('y (m)');
title('original dataset');
subplot(1,2,2)
mesh(X,Y,Z); view(2); grid on;
title('interpolated grid');
axis equal tight; xlabel('x (m)'); yticklabels('');

Additional code, using alphaShape

shp = alphaShape(x,y,dx);
figure
plot(shp)
title("Alpha Shape, alpha=5")

Use the alpha shape to mask values in Z.

Z(~inShape(shp,X,Y)) = NaN;
mesh(X,Y,Z); view(2); grid on;
title('interpolated grid');
axis equal tight; xlabel('x (m)'); yticklabels('');

[Lorenzo Melito]

Thank you very much for pointing me at a very useful function I didn't know of, Mr. Eddins! alphaShape works great for my application. Thanks again and kind regards!