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],.开发者_JAVA百科...]
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)?

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

This doesn't sound so much to me like interpolation as finding which trianglular face the point is interior to. Check this site out for a way to test each triangular face. Your function will just determine which face it is inside and return the corresponding value. Of course, if you have a lot of faces or if you are doing this a whole lot, then you will want to find ways to optimize this (store the farthest + and - points for each triangle in both the x and y directions and store them with the face for instance. If the point isn't within this bounding box, then you may as well not check if it is interior to the triangle).

I really doubt that you will find something that is built-in to Matlab or scipy to do just what you want, but I could be wrong.

you might want to use the bridge CGAL-python module; If I remember correctly, CGAL provides functionality for triangle search. It probably works most efficiently, however, if you construct the triangulation incrementally using their built-in representation. For a quick-and-dirty one off, you could just find the closest mesh vertex to the query point, either by a Voronoi diagram (the functionality for this in Matlab is not great), or, for a single query, computing all the distances, and finding the minimum, then search all triangles with that vertex.

Matlab has the built-in function inpolygon, which allows you to test whether you're inside a triangle. I do not know of a function that would identify inside which face you are.

If you were to write such a function, I'd first test which vertex your point is closest to, and then evaluate inpolygon for all faces that share the vertex until you find the match. This should be reasonably fast.

Take a look at matplotlib.delaunay.interpolate, which is a well-documented wrapper for C code.
(However class LinearInterpolator there says "At the moment, only regular rectangular grids are supported for interpolation.")