matlab - 在面上插值分段常数的二维数据

Question

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

我正在寻找一种方法来根据这些数据构造函数 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)中实现了这个功能？

Answer 1

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

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

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

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

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

• 输入检查:代码目前假设faces是一个N-by-3矩阵，vertices是一个M-by-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