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 < dmin
         dmin = d;
      end
   end
   if dmin >= 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 full mesh, and here the same mesh after the undesired points are removed. In both pictures the red dots represent the initial data set the mesh is interpolated from.

My main concern here is that while this code works fine with a limited number of points in the starting data set, it gets incredibly cumbersome when the number of data points is in the order of hundreds of thousands, which is sadly my case. I can let it be if necessary, but I feel like I should make the code a bit more efficient, since I may have to run it many times. Do you have any suggestion on how to improve its efficiency? Any input will be greatly appreciated.

I am using MATLAB 2017a.