在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)
triVertex = vertices(faces(:,2),:);
triLegLeft = vertices(faces(:,1),:)-triVertex;
triLegRight = vertices(faces(:,3),:)-triVertex;
C1 = sum(triLegLeft.*triLegRight,2);
C2 = sum(triLegLeft.^2,2);
C3 = sum(triLegRight.^2,2);
triBoundary = max(C2,C3);
scale = C1.^2-C2.*C3;
C1 = C1./scale;
C2 = C2./scale;
C3 = C3./scale;
f = @evaluate_function;
function val = evaluate_function(x,y)
w = [x-triVertex(:,1) y-triVertex(:,2)];
nearIndex = find(sum(w.^2,2) <= triBoundary);
if isempty(nearIndex)
val = NaN;
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;
else
val = values(nearIndex(inIndex));
end
end
end