poster
2019-12-10, 20:30
我有一个由两个变量描述的不规则网格:一个用于存储构成每个面的顶点索引的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的建议 (https://stackoverflow.com/questions/2300103/interpolating-2d-data-that-is-piecewise-constant-on-faces/2300736#2300736)使用INPOLYGON (http://www.mathworks.com/access/helpdesk/help/techdoc/ref/inpolygon.html)函数构建自己的算法,但是您可以使用一些标准算法来查找点是否在多边形内,从而自己创建一个更快的实现。
前一阵子,我写了自己的代码来查找线段和一组3-D中的三角形表面之间的交点,并且我发现此softsurfer链接 (http://softsurfer.com/Archive/algorithm_0105/algorithm_0105.htm)对于实现算法最有用。我的案子比你的案子复杂。由于您使用的是2-D,因此可以忽略链接的第一部分,即查找线段与三角形平面相交的点。
我在下面提供了我的MATLAB代码的简化版本,供您使用。函数interpolate会将您的faces , vertices和values矩阵作为输入,并返回可以在给定(x,y)点求值的函数句柄f ,以获取边界三角形内的分段值。以下是此代码的一些功能:
评估f时将执行的处理包含在嵌套函数 (http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_prog/f4-39683.html) evaluate_function 。该功能可以访问在其它变量interpolate ,所以一些所需的三角测试变量预先计算以便evaluate_function运行尽可能快。
如果您有很多三角形,则测试点是否在所有三角形内可能会很昂贵。因此,代码首先查找在您的点的给定半径(即最长的三角形的长度)内的三角形。仅测试这些附近的三角形,以查看该点是否在其中。
如果点未落在任何三角形区域内,则f返回NaN (http://www.mathworks.com/access/helpdesk/help/techdoc/ref/nan.html)的值。
代码中未包含某些内容,您可能希望根据其用途将这些内容添加进去:
输入检查:代码当前假设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
我正在寻找一种从该数据构造函数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的建议 (https://stackoverflow.com/questions/2300103/interpolating-2d-data-that-is-piecewise-constant-on-faces/2300736#2300736)使用INPOLYGON (http://www.mathworks.com/access/helpdesk/help/techdoc/ref/inpolygon.html)函数构建自己的算法,但是您可以使用一些标准算法来查找点是否在多边形内,从而自己创建一个更快的实现。
前一阵子,我写了自己的代码来查找线段和一组3-D中的三角形表面之间的交点,并且我发现此softsurfer链接 (http://softsurfer.com/Archive/algorithm_0105/algorithm_0105.htm)对于实现算法最有用。我的案子比你的案子复杂。由于您使用的是2-D,因此可以忽略链接的第一部分,即查找线段与三角形平面相交的点。
我在下面提供了我的MATLAB代码的简化版本,供您使用。函数interpolate会将您的faces , vertices和values矩阵作为输入,并返回可以在给定(x,y)点求值的函数句柄f ,以获取边界三角形内的分段值。以下是此代码的一些功能:
评估f时将执行的处理包含在嵌套函数 (http://www.mathworks.com/access/helpdesk/help/techdoc/matlab_prog/f4-39683.html) evaluate_function 。该功能可以访问在其它变量interpolate ,所以一些所需的三角测试变量预先计算以便evaluate_function运行尽可能快。
如果您有很多三角形,则测试点是否在所有三角形内可能会很昂贵。因此,代码首先查找在您的点的给定半径(即最长的三角形的长度)内的三角形。仅测试这些附近的三角形,以查看该点是否在其中。
如果点未落在任何三角形区域内,则f返回NaN (http://www.mathworks.com/access/helpdesk/help/techdoc/ref/nan.html)的值。
代码中未包含某些内容,您可能希望根据其用途将这些内容添加进去:
输入检查:代码当前假设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