Removing undesired points from a mesh

I have a 3D data set made of a certain number of points (x,y,z) that cover a certain region of space. Using the scatteredInterpolant object I can interpolate this data set over a grid to produce a rectangular mesh. Note that the mesh may extend to regions that are not defined by the data set; in fact, after the mesh generation I need to remove the part of the mesh that is extrapolated away from the data set (replacing its values with NaN, for example) in order to retain only the mesh generated between the data points.

I came up with the following MATLAB script to solve this problem in a very naive way. Considering a single mesh point (xq,yq), I evaluate the minimum distance between this point and the data set; if this distance is greater than a certain threshold, then the corresponding interpolated value (zq) is set to NaN.

%% Data set (x,y,z)
x = [3 3 3 4 4 4 4 4 5 5 5 5 5]';
y = [1 2 3 0 1 2 3 4 0 1 2 3 4]';
z = [.5 .505 .51 .51 .51 .51 .51 .515 .535 .528 .53 .53 .53]';

%% Interpolant
F = scatteredInterpolant(x,y,z,'natural');

%% Mesh generation (xq,yq,zq)
delta = 0.5;
ti = 0:delta:5;
si = 0:delta:4;
[xq,yq] = meshgrid(ti,si);
zq = F(xq,yq);

%% Replacing undesired values with NaN
thresh = 1;
n = length(ti) * length(si);
m = length(x);
xqcol = reshape(xq,[n,1]);
yqcol = reshape(yq,[n,1]);
zqcol = reshape(zq,[n,1]);
tab = [xqcol yqcol zqcol];

for i = 1:n
   dmin = 10^32;
   for k = 1:m
      diffx = tab(i,1) - x(k);
      diffy = tab(i,2) - y(k);
      d = sqrt(diffx^2 + diffy^2);
      if d = thresh
      tab(i,3) = NaN;
   end
end

zqwork = tab(:,3);
zq2 = reshape(zqwork,[size(zq,1),size(zq,2)]);

%% Plotting
figure
plot3(x,y,z,'.r','MarkerSize',10);   % data set
grid on; axis([0 5 0 4 0.46 0.54]);
hold on;
% mesh(xq,yq,zq); view(3);           % full mesh
mesh(xq,yq,zq2); view(3);            % mesh with undesired points removed

This code gets the job done (to my utter surprise, since I am not really proficient with Matlab!). Here you can find a fu

MATLAB's documentation about the scatterInterpolant function actually mentions in passing how one would do this in https://www.mathworks.com/help/matlab/ref/scatteredinterpolant-object.html#bvkm1tv-4. In short, you would set F.ExtrapolationMethod to None so that F returns NaN for points that do not lie inside the convex hull of the points defined by x and y. This way, you can avoid writing for loops and a bunch of reshape commands, which can take a while for large vectors and matrices.

x = [3 3 3 4 4 4 4 4 5 5 5 5 5]';
y = [1 2 3 0 1 2 3 4 0 1 2 3 4]';
z = [.5 .505 .51 .51 .51 .51 .51 .515 .535 .528 .53 .53 .53]';

F = scatteredInterpolant(x,y,z,'natural');

delta = 0.5;
ti = 0:delta:5;
si = 0:delta:4;
[xq,yq] = meshgrid(ti,si);
zq = F(xq,yq);

figure
hold on
mesh(xq,yq,zq)
plot3(x,y,z,'r.','MarkerSize',10)
hold off
view(3)
axis([0 5 0 4 0.46 0.54])

F.ExtrapolationMethod = 'none';
zqe = F(xq,yq);

figure
hold on
mesh(xq,yq,zqe)
plot3(x,y,z,'r.','MarkerSize',10)
hold off
view(3)
axis([0 5 0 4 0.46 0.54])

Plot of all points, including extrapolated ones:

Plot of points inside convex hull, excluding extrapolated ones:

StackExchange Code Review Q#160735, answer score: 2