logo标签列表

在面上分段常数的情况下插值2D数据

matlabscipyinterpolationpiecewise
3

我有一个不规则的网格，由两个变量描述 - 存储构成每个面的顶点索引的面数组和存储每个顶点坐标的 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）中实现这个函数的方法？
5个回答
1

在MATLAB中没有内置的函数可以实现您想要的功能。您可以使用INPOLYGON函数，如Jonas所建议的，构建自己的算法，但是您也可以使用一些标准算法来创建更快的实现，以查找点是否在多边形内。

一段时间以前，我编写了自己的代码，用于在3D中找到线段和一组三角形表面之间的交点，并且我发现this softsurfer link对于实现该算法非常有帮助。我的情况比您的更复杂。由于您正在使用2D，因此可以忽略链接的第一部分，该部分介绍了查找线段与三角形平面相交的点的方法。

我已经为您提供了一个简化版本的MATLAB代码。函数interpolate将以您的faces、vertices和values矩阵作为输入，并返回一个函数句柄f，该句柄可以在给定的（x，y）点处进行评估，以获取边界三角形内的分段值。以下是此代码的一些特点：
  • 当你评估f时将执行的处理包含在嵌套函数evaluate_function中。该函数可以访问interpolate中的其他变量，因此为了使evaluate_function运行尽可能快，需要预先计算一些用于三角形内测试所需的变量。
  • 在有大量三角形的情况下，测试您的点是否在其中所有三角形内可能会很昂贵。因此，代码首先找到距离您的点给定半径（即三角形最长边的长度）内的三角形。只有这些附近的三角形才会被测试以查看点是否在其中。
  • 如果一个点不在任何三角形区域内，则f返回NaN的值。

代码中没有包含一些内容，您可能想根据使用情况添加：

  • 输入检查： 该代码当前假定faces是一个N x 3矩阵，vertices是一个M x 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
1
这对我来说不像是插值，而更像是找到点所在的三角形面。可以查看this site以了解测试每个三角形面的方法。您的函数将确定它在哪个面内并返回相应的值。当然，如果您有很多面或者您经常这样做，那么您需要找到优化的方法（例如，在x和y方向上存储每个三角形的最远+和-点，并将它们与面一起存储。如果该点不在此边界框内，则您可能不需要检查它是否在三角形内）。
我真的怀疑您能够在Matlab或scipy中找到内置的功能来完成您想要的事情，但我可能错了。
1

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

1
Matlab有内置函数inpolygon，可让您测试是否在三角形内。我不知道有哪个函数可以确定您在哪个面内。
如果您要编写这样的函数，我将首先测试您的点最靠近哪个顶点，然后对共享该顶点的所有面进行inpolygon评估，直到找到匹配项。这应该是相当快的。
1

看一下 matplotlib.delaunay.interpolate，这是一个有良好文档的 C 代码包装器。
（然而那里的 class LinearInterpolator 说：“目前只支持规则矩形网格进行插值。”）

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接