插补在表面上分段恒定的2d数据

[poster]

我有一个由两个变量描述的不规则网格：一个用于存储构成每个面的顶点索引的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）中实现？



回答：

MATLAB中没有内置函数可以执行您想要的操作。您可以按照Jonas的建议使用INPOLYGON函数构建自己的算法，但是您可以使用一些标准算法来查找点是否在多边形内，从而自己创建一个更快的实现。

前一阵子，我写了自己的代码来查找线段和一组3-D中的三角形表面之间的交点，并且我发现此softsurfer链接对于实现算法最有用。我的案子比你的案子复杂。由于您使用的是2-D，因此可以忽略链接的第一部分，即查找线段与三角形平面相交的点。

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

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

代码中未包含某些内容，您可能希望根据其用途将这些内容添加进去：

• 输入检查：代码当前假设faces是一个N×3矩阵， vertices是一个M×2矩阵，而values是一个长度为N的向量。您可能需要添加错误检查，以确保输入符合这些要求，并抛出错误以指示何时其中一个或多个不正确。
• 退化三角形检查：由您的faces和vertices输入定义的一个或多个三角形可能会退化（即，它们的面积可能为0）。当两个或更多个三角形顶点是相同的精确点时，或者当三角形的所有三个顶点都位于一条直线上时，就会发生这种情况。您可能想添加一个检查，在计算f时将忽略这些三角形。
• 处理边缘情况：一个点可能会终止于两个或多个三角形区域的边缘。因此，您必须决定该点将采用的值（即，最大的面值，面值的平均值等）。对于像这样的边缘情况，下面的代码将自动选择接近于您的faces变量中的faces列表开头的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) = 0) & (t >= 0) & (s+t