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

14 views (last 30 days)
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?

Accepted Answer

Steve Eddins
Steve Eddins on 26 Jan 2022
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('');
  1 Comment
Lorenzo Melito
Lorenzo Melito on 26 Jan 2022
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!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!