对面上分段常数的二维数据进行插值

发布于 2024-08-22 09:51:14 字数 554 浏览 12 评论 0原文

我有一个由两个变量描述的不规则网格 - 一个 faces 数组，用于存储构成每个面的顶点的索引，以及一个 verts 数组，用于存储每个顶点的坐标。我还有一个函数，假设每个面上都是分段常数，并且它以每个面值数组的形式存储。

我正在寻找一种从这些数据构造函数 f 的方法。大致如下：

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]

计算 f(x,y) 的手动方法是找到点 x,y 所在的相应面，并返回存储在该面上的值。是否有一个函数已经在 scipy （或 matlab ）中实现了这个功能？

原文

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)?

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

猛虎独行 2024-08-29 09:51:14

MATLAB 中没有内置函数可以完成您想要的操作。您可以使用函数 INPOLYGON 正如乔纳斯建议的，但是您可以使用一些标准算法自行创建更快的实现来查找点是否在多边形内。

不久前，我编写了自己的代码，用于在 3D 中查找线段和一组三角形表面之间的交点，并且我发现此 softsurfer 链接对于实现该算法最有帮助。我的情况比你的复杂。由于您正在处理二维，因此您可以忽略有关查找线段与三角形平面相交点的链接的第一部分。

我在下面提供了 MATLAB 代码的简化版本供您使用。函数interpolate 会将您的面、顶点 和值 矩阵作为输入并返回函数句柄f 可以在给定的 (x,y) 点求值以获得边界三角形内的分段值。以下是此代码的一些功能：

当您计算 f 时将完成的处理包含在嵌套函数 evaluate_function。该函数可以访问 interpolate 中的其他变量，因此预先计算了三角形内测试所需的许多变量，以便 evaluate_function 尽可能快地运行。
如果您有很多三角形，测试您的点是否在所有三角形内可能会很昂贵。因此，代码首先查找点的给定半径（即三角形最长边的长度）内的三角形。仅测试这些附近的三角形以查看该点是否在其中。
如果点未落在任何三角形区域内，则值为 NaN 由 f 返回。

有些内容未包含在代码中，您可能需要根据您的用途添加这些内容：

输入检查：代码当前假设faces 是一个 N×3 矩阵，vertices 是一个 M×2 矩阵，values 是一个长度为 N 的向量。您可能希望添加错误检查以确保输入符合这些要求，并在其中一个或多个输入不正确时抛出错误。
简并三角形检查：由您的面和顶点输入定义的一个或多个三角形可能是简并的（即它们可能具有面积为0）。当三角形的两个或多个顶点精确地位于同一点时，或者当三角形的所有三个顶点都位于一条直线上时，就会发生这种情况。您可能希望添加一个检查，在评估 f 时忽略此类三角形。
处理边缘情况：一个点可能最终位于两个或多个三角形区域的边缘上。因此，您必须决定该点将采用什么值（即面值中的最大值、面值的平均值等）。对于像这样的边缘情况，下面的代码将自动选择更接近 faces 变量中面部列表开头的面部值。

最后，这是代码：

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

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

回复收藏 0 原文

能怎样 2024-08-29 09:51:14

对我来说，这听起来不像插值，而是找到该点位于哪个三角形面的内部。查看此网站，了解测试每个三角形面的方法。您的函数将仅确定它位于哪个面内并返回相应的值。当然，如果你有很多面或者你做了很多次，那么你会想要找到优化它的方法（存储每个三角形在 x 和 y 方向上最远的 + 和 - 点，并存储例如，如果该点不在该边界框内，那么您也可以不检查它是否在三角形内部）。

我真的怀疑你是否会找到 Matlab 或 scipy 内置的东西来完成你想要的事情，但我可能是错的。

回复收藏 0 原文

风吹雪碎 2024-08-29 09:51:14

您可能想使用桥 CGAL-python 模块；如果我没记错的话，CGAL 提供了三角形搜索的功能。但是，如果您使用其内置表示逐步构建三角剖分，它可能会最有效地工作。对于快速而肮脏的一件事，您可以通过 Voronoi 图（Matlab 中的功能不是很好）找到距离查询点最近的网格顶点，或者对于单个查询，计算所有距离，找到最小值，然后搜索具有该顶点的所有三角形。

回复收藏 0 原文