Interpolating 2d data that is piecewise constant on faces

I have an irregular mesh which is described by two variables - a faces array that stores the indices of the vertices that constitute each face, and a verts array that stores the coordinates of each vertex. I also have a function that is assumed to be piecewise constant over each face, and it is stored in the form of an array of values per face.

I am looking for a way to construct a function f from this data. Something along the following lines:

faces = [[0,1,2], [1,2,3], [2,3,4] ...]
verts = [[0,0], [0,1], [1,0], [1,1],....]
vals = [0.0, 1.0, 0.5, 3.0,....]

f = interpolate(faces, verts, vals)

f(0.2, 0.2) = 0.0 # point inside face [0,1,2]
f(0.6, 0.6) = 1.0 # point inside face [1,2,3]

The manual way of evaluating f(x,y) would be to find the corresponding face that the point x,y lies in, and return the value that is stored in that face. Is there a function that already implements this in scipy (or in matlab)?

Solution

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