8

假设我有 n 个点定义 z 轴上的表面

f(x1,y1) = 10
f(x2,y2) = 12
f(x3,y3) = 5
f(x4,y4) = 2
...
f(xn,yn) = 21

现在我希望能够近似 f(x,y)。我正在寻找一种线性算法,尤其是样条近似算法。一个示例算法或至少一些指针会很棒。

4

3 回答 3

4

这是对进行线性近似的方法的模糊描述。

  1. 确定您的点的Voronoi 图(对于平面中的每个点,找到最近的(x_i,y_i)
  2. 取它的对偶以获得Delaunay 三角剖分:连接(x_i,y_i)(x_j,y_j)如果有一条线段的点,那么(x_i,y_i)(x_j,y_j)是等距的(并且比任何其他对更近)。
  3. 在每个三角形上,找到通过三个顶点的平面。这是您需要的线性插值。

下面实现了 Python 中的前两个步骤。网格的规律性可能会让你加快速度(它也可能会打乱三角测量)。

import itertools

""" Based on http://stackoverflow.com/a/1165943/2336725 """
def is_ccw(tri):
    return ( ( tri[1][0]-tri[0][0] )*( tri[1][1]+tri[0][1] )
            + ( tri[2][0]-tri[1][0] )*( tri[2][1]+tri[1][1] )
            + ( tri[0][0]-tri[2][0] )*( tri[0][1]+tri[2][1] ) ) < 0

def circumcircle_contains_point(triangle,point):
    import numpy as np
    matrix = np.array( [ [p[0],p[1],p[0]**2+p[1]**2,1] for p in triangle+point ] )
    return is_ccw(triangle) == ( np.linalg.det(matrix) > 0 )

triangulation = set()
"""
A triangle is in the Delaunay triangulation if and only if its circumscribing
circle contains no other point.  This implementation is O(n^4).  Faster methods
are at http://en.wikipedia.org/wiki/Delaunay_triangulation
"""
for triangle in itertools.combinations(points,3):
    triangle_empty = True
    for point in points:
        if point in triangle:
            next
        if circumcircle_contains_point(triangle,point):
            triangle_empty = False
            break
    if triangle_empty:
        triangulation.add(triangle)
于 2014-05-05T21:43:37.930 回答
2

对不规则 2D 数据进行插值并不容易。我知道没有真正的样条泛化到不规则的二维。

除了基于三角测量的方法之外,您还可以查看 Barnes ( http://en.wikipedia.org/wiki/Barnes_interpolation ) 和反距离加权 ( http://en.wikipedia.org/wiki/Inverse_distance_weighting ),或者更一般地说,RBF(http://en.wikipedia.org/wiki/Radial_basis_functions)。

如果您的点非常不均匀地分布(密集簇),则可能有必要使函数的大小自适应,或者求助于近似而不是插值。

于 2014-05-06T08:26:24.620 回答
1

您可以将您的点用作 Bezier(或 Bspline)曲面的控制点,尤其是在对平面(xi, yi)中的矩形进行采样时。XY在这方面,不涉及任何拟合。

您将获得的表面将位于点的凸包中,并将与 边界处的点相交(插值){xi, yi}

如果您想尝试一下,这个论坛帖子似乎包含简单的代码,如果您没有Matlab,您可以使用GuIRITMatlab来做同样的事情(尽管它需要弄清楚程序的文件格式)。

于 2012-03-11T17:36:50.060 回答