Removing undesired points from a mesh

Question

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.

Answer 1

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:

enter image description here

Plot of points inside convex hull, excluding extrapolated ones:

enter image description here

• 19lorenz88: thank you very much for your input! I totally missed that option while I was reading the documentation. I should have read it more thoroughly. But now I realize my example was a bit partial. It seems like the extrapolation can be avoided only outside of the convex hull of the data: what if the starting data has any "holes" inside the hull?
• edwinksl: A convex hull can't have holes though, so I am not sure what you mean.
• 19lorenz88: I hope this gif clarifies what I mean. Suppose the data set has some sort of "hole", a portion of space in which there is no point (in the example three red points are missing at the center of the data set). scatteredInterpolant generates an interpolating mesh over the whole domain [xq yq], but the option ExtrapolationMethod set to none excludes only the points outside the convex hull, while the mesh inside the "hole" is retained.
• 19lorenz88: I can use my script to get rid of the mesh inside the hole (I tested it and it works pretty fine for my purposes), but when it comes to processing a data set of thousands of points the task takes quite a while. I've skimmed through the documentation and I could not find any ready option for scatteredInterpolant to do the job. Maybe another function will do?
• edwinksl: Even if you remove those 3 points, the convex hull remains the same as before. Since any point in the convex hull is considered to be an interpolation while any point outside it is considered an extrapolation, the situation is basically the same as before (but not quantitatively the same because removing these 3 points does affect how the interpolant behaves in their neighborhood).