There's no built-in function in MATLAB that will do what you want. You could likely build your own algorithm using the function INPOLYGON as suggested by Jonas, but you may be able to create a faster implementation yourself using some standard algorithms for finding whether a point is within a polygon.
A while back I wrote my own code for finding the intersection points between a line segment and a set of triangular surfaces in 3-D, and I found this softsurfer link to be the most helpful for implementing the algorithm. My case was more complicated than yours. Since you are working in 2-D, you can ignore the first section of the link about finding the point where the segment intersects the plane of the triangle.
I've included a simplified version of my MATLAB code below for you to use. The function interpolate will take your faces, vertices, and values matrices as inputs and return a function handle f that can be evaluated at a given (x,y) point to get the piecewise value within the bounding triangle. Here are a few features of this code:
- The processing that will be done when you evaluate
f is contained in the nested function evaluate_function. This function has access to the other variables in interpolate, so a number of variables needed for the in-triangle test are precomputed so that evaluate_function runs as fast as possible.
- In cases where you have a lot of triangles, testing if your point is inside all of them could be expensive. Therefore, the code first finds triangles that are within a given radius (i.e. the length of the triangles longest leg) of your point. Only these nearby triangles are tested to see if the point is within them.
- If a point does not fall inside any triangular region, a value of NaN is returned by
f.
There are some things that are not included in the code, which you may want to add in depending on what you are using it for:
- Input check: The code currently assumes that
faces is an N-by-3 matrix, vertices is an M-by-2 matrix, and values is a length N vector. You would likely want to add error checking to make sure that the inputs conform to these requirements, and throw an error indicating when one or more of them is incorrect.
- Degenerate triangle check: It's possible that one or more of the triangles defined by your
faces and vertices inputs may be degenerate (i.e. they may have an area of 0). This occurs when two or more of the triangles vertices are the same exact point, or when all three vertices of the triangle lie in a straight line. You would likely want to add a check that will ignore such triangles when it comes to evaluating f.
- Handling edge cases: It's possible that a point could end up on the edge of two or more triangular regions. You would therefore have to decide what value the point will take on (i.e. the largest of the face values, the mean of the face values, etc.). For edge cases like this, the code below will automatically pick the value of the face that is closer to the beginning of the list of faces in your
faces variable.
Finally, here's the code:
function f = interpolate(faces,vertices,values)
%# Precompute some data (helps increase speed):
triVertex = vertices(faces(:,2),:); %# Triangles main vertices
triLegLeft = vertices(faces(:,1),:)-triVertex; %# Triangles left legs
triLegRight = vertices(faces(:,3),:)-triVertex; %# Triangles right legs
C1 = sum(triLegLeft.*triLegRight,2); %# Dot product of legs
C2 = sum(triLegLeft.^2,2); %# Squared length of left leg
C3 = sum(triLegRight.^2,2); %# Squared length of right leg
triBoundary = max(C2,C3); %# Squared radius of triangle boundary
scale = C1.^2-C2.*C3;
C1 = C1./scale;
C2 = C2./scale;
C3 = C3./scale;
%# Return a function handle to the nested function:
f = @evaluate_function;
%# The nested evaluation function:
function val = evaluate_function(x,y)
w = [x-triVertex(:,1) y-triVertex(:,2)];
nearIndex = find(sum(w.^2,2) <= triBoundary); %# Find nearby triangles
if isempty(nearIndex)
val = NaN; %# Return NaN if no triangles are nearby
return
end
w = w(nearIndex,:);
wdotu = sum(w.*triLegLeft(nearIndex,:),2);
wdotv = sum(w.*triLegRight(nearIndex,:),2);
C = C1(nearIndex);
s = C.*wdotv-C3(nearIndex).*wdotu;
t = C.*wdotu-C2(nearIndex).*wdotv;
inIndex = find((s >= 0) & (t >= 0) & (s+t <= 1),1);
if isempty(inIndex)
val = NaN; %# Return NaN if point is outside all triangles
else
val = values(nearIndex(inIndex));
end
end
end