基于移动最小二乘法的曲线曲面拟合(python语言实现)

1.移动最小二乘法

上篇论文采用最小二乘法来拟合曲线，如果离散数据量比较大，形状复杂，还需要分段拟合和平滑化，因此采用移动最小二乘法进行曲线拟合，可以克服上面的缺点，还具有一些优点；
移动最小二乘法与传统的最小二乘法相比，有两个比较大的改进：

（ 1）拟合函数的建立不同。这种方法建立拟合函数不是采用传统的多项式或其它函数，而是由一个系数向量 a(x)和基函数 p(x)构成， 这里 a(x)不是常数，而是坐标 x 的函数。
（ 2）引入紧支（ Compact Support）概念，认为点 x 处的值 y 只受 x
附近子域内节点影响，这个子域称作点 x 的影响区域， 影响区域外的节点对 x的取值没有影响。在影响区域上定义一个权函数w(x)， 如果权函数在整个区域取为常数， 就得到传统的最小二乘法。

参考自《基于移动最小二乘法的曲线曲面拟合-曾清红》

2.拟合函数的建立

在拟合区域的一个局部子域上， 拟合函数 f (x)表示为：

式中 为待求系数，它是坐标x的函数。称为基函数。它是一个k阶完备的多项式，m是基函数的项数，

对于一维问题 ：
基函数可以为 p（x）=$[1,x,x^2,,,,x^m]$

二维问题可以为： 线性基 p（x）=$[1.x.y]^T$, m=3 二次基 p(x) =$[1,x,y,x^2,xy,y^2]^T$ m=6

这是为在阅读文献时的疑惑，因为我解决的是一维问题，所以不需要二维的基函数。
在移动最小二乘近似中, 系数 $a_i(x)$ 是通过令近似函数 u(x) 在点 x 的邻域 内各节点误差的加权平方和为最小来确定的


式中 n 为点 x 的邻域 内所包含的节点数.,$w(x) = (x -x_I)$ 称为节点 $x_I$ 处的权函数, 它在节点 xI 周围的一个有限区域中大于零, 而在该区域外为零 . 权函数的定义表明, 只有在节点 xI 的影响域范围内的节点才对该点的近似函数产生影响.

这里对支撑域进行说明：
如图：

将整个x范围划分为若干个区域，每个区域包含若干个x点，那么并且规定其中一点为标准点，其他点为参考点。
参考点与标准点的距离作为权函数的参数。得出权重。

3.权函数

权函数在移动最小二乘法中起着非常重要的作用。移动最小二乘法中的权函数 $w(x-x_I)$应该具有紧支性，也就是权函数在 x
的一个子域内不等于零， 在这个子域之外全为零， 这个子域称为权函数的支持域(即 x 的影响区域)。一般选择圆形作为权函数的支持域(见图其半径记为 $s_{max}$。 由于权函数的紧支性，只有这些包含在影响区域内的数据点对点 x 的取值有影响权函数 $w(x-x_I)$应该是非负的，并且随着$||x - x_i||^2$ 的增加单调递减。权函数还应具有一定的光滑性，因为拟合函数会继承权函数的连续性：如果权函数$w(x-x_I)$是 C1 阶连续的，则拟合函数也是 C1 阶连续的。常用的权函数是样条函数

4

3 法方程的推导
对于任意函数 h(x) 和 g(x), 引入记号：

那么：
公式4可以写为：

写成矩阵形式：

由上面的法方程, 解得 a(x).
然后求解得出A(x),B(x) 求解得出$\alpha$(x)

4.拟合流程

这里说明一下为什么要网格化，网格化主要是选取标准点，并以标准点来划分支撑域，确定支撑域半径和支撑域内的节点
x。
我仍然以上篇最小二乘法的数据点为例，通过代码编写移动最小二乘法的方法：

#主题部分
X=np.arange(-0.9,0.9,0.05)
# 数据点x个数
M=len(x)
# 基函数个数
N=2
p=np.zeros((M,2))
Y=[]
for XX in X:
    w = np.zeros((M,1))
    d=0.1 # 影响区域的半径
    for i in range(0,M):
        w[i]=W_fun(d, x[i], XX)
        p[i][0]=1
        p[i][1]=x[i]
    A=fun_A(x,w,p)
    B=fun_B(y,w,p)
    a=np.linalg.solve(A,B)
    Y.append(a[0]+a[1]*XX)


----------


#其他函数部分
# 权函数
def W_fun(d,x,X):
    s=abs(x-X)/d
    if (s<=0.5):
        return (2/3)-4*s**2+4*s**3
    elif(s<=1):
        return (4/3)-4*s+4*s**2-(4/3)*s**3
    else:
        return 0
# 权函数记号(pm,pn)的计算
def pm_pn(w,x,p,m,n):
    # x为数据点，w为支撑域的权重，M为数据点个数 p1,p2为传入的数值
    pmn=0
    M=len(x)
    # i代表数据点,m n代表(pm,pn)的下标
    for i in range(M):
        pmn=pmn+w[i]*p[i][m]*p[i][n]
    return float(pmn)
# B矩阵的建立
def fun_B(u,w,p):
    pumI=0
    M=len(u) #数据点个数
    m=p.shape[1] # 基函数个数
    B=[]
    for j in range(m):
        for i in range(M):
            pumI=pumI+w[i]*p[i][j]*u[i]
        B.append(float(pumI))
    return B
 # A矩阵的建立
def fun_A(x,w,p):
    M=len(x)函数
    m=p.shape[1]
    A=[]
    for mm in range(m):
        matA=[]
        for nn in range(m):
            pmn=pm_pn(w,x,p,mm,nn)
            matA.append(pmn)
        A.append(matA)
    return A

结果

5.结果

绿色为移动最小二乘法，红色为最小二乘法。

参考文献：
1基于移动最小二乘法的曲线曲面拟合-曾清红
2移动最小二乘法在多功能传感器数据重构中的应用-刘丹
3 移动最小二乘法（MLS）曲线曲面拟合C++代码实现
https://blog.csdn.net/liumangmao1314/article/details/54179526