3

我有一个由两个变量描述的不规则网格 - 一个 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 ）中实现此功能的函数？

4

5 回答 5

1

MATLAB 中没有内置函数可以满足您的需求。您可能会按照 Jonas的建议使用函数INPOLYGON构建自己的算法，但您可以使用一些标准算法自己创建更快的实现来查找点是否在多边形内。

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

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

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

代码中未包含一些内容，您可能需要添加这些内容，具体取决于您使用它的目的：

  • 输入检查：代码当前假设faces是一个 N×3 矩阵，vertices是一个 M×2 矩阵，并且values是一个长度为 N 的向量。您可能希望添加错误检查以确保输入符合这些要求，并在其中一个或多个不正确时引发错误。
  • 退化三角形检查：由您faces和vertices输入定义的一个或多个三角形可能是退化的（即它们的面积可能为 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
于 2010-03-08T16:17:50.493 回答
1

看看matplotlib.delaunay.interpolate，它是一个有据可查的 C 代码包装器。
（但是class LinearInterpolator那里说“目前，仅支持常规矩形网格进行插值。”）

于 2010-03-07T20:08:49.927 回答
1

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

于 2010-02-20T02:01:02.353 回答
1

Matlab 有内置函数inpolygon，它允许你测试你是否在一个三角形内。我不知道有什么功能可以识别你在哪张脸里面。

如果您要编写这样的函数，我将首先测试您的点最接近哪个顶点，然后评估共享该顶点的所有面的多边形，直到找到匹配项。这应该相当快。

于 2010-02-20T02:01:14.533 回答
1

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

我真的怀疑你会找到内置于 Matlab 或 scipy 的东西来做你想做的事，但我可能是错的。

于 2010-02-19T23:22:19.880 回答