2.移動最小二乘法
轉(zhuǎn)載自：基于移動最小二乘法的曲線曲面擬合--一碗竹葉青
如果離散數(shù)據(jù)量比較大，形狀復(fù)雜，用傳統(tǒng)最小二乘法會很奇怪。
所以又有了改良型的基于移動最小二乘法。
區(qū)別：
（ 1）擬合函數(shù)的建立不同。這種方法建立擬合函數(shù)不是采用傳統(tǒng)的多項(xiàng)式或其它函數(shù)，而是由一個系數(shù)向量 a(x)和基函數(shù) p(x)構(gòu)成， 這里 a(x)不是常數(shù)，而是坐標(biāo) x 的函數(shù)。
（ 2）引入緊支（ Compact Support）概念，認(rèn)為點(diǎn) x 處的值 y 只受 x
附近子域內(nèi)節(jié)點(diǎn)影響，這個子域稱作點(diǎn) x 的影響區(qū)域， 影響區(qū)域外的節(jié)點(diǎn)對 x的取值沒有影響。在影響區(qū)域上定義一個權(quán)函數(shù)w(x)， 如果權(quán)函數(shù)在整個區(qū)域取為常數(shù)， 就得到傳統(tǒng)的最小二乘法。參考自《基于移動最小二乘法的曲線曲面擬合-曾清紅》
1.擬合函數(shù)的建立
在擬合區(qū)域的一個局部子域上， 擬合函數(shù) f (x)表示為：
式中
這個基函數(shù)我還真的看不懂，但大家都是這樣設(shè)的
在移動最小二乘近似中, 系數(shù) ai(x) 是通過令近似函數(shù) u(x) 在點(diǎn) x 的鄰域內(nèi)各節(jié)點(diǎn)誤差的加權(quán)平方和為最小來確定的
式中 n 為點(diǎn) x 的鄰域 內(nèi)所包含的節(jié)點(diǎn)數(shù).,w(x)=(x?xI) 稱為節(jié)點(diǎn) xI 處的權(quán)函數(shù), 它在節(jié)點(diǎn) xI 周圍的一個有限區(qū)域中大于零, 而在該區(qū)域外為零 . 權(quán)函數(shù)的定義表明, 只有在節(jié)點(diǎn) xI 的影響域范圍內(nèi)的節(jié)點(diǎn)才對該點(diǎn)的近似函數(shù)產(chǎn)生影響.
2.支撐域：
將整個x范圍劃分為若干個區(qū)域，每個區(qū)域包含若干個x點(diǎn)，那么并且規(guī)定其中一點(diǎn)為標(biāo)準(zhǔn)點(diǎn)，其他點(diǎn)為參考點(diǎn)。
參考點(diǎn)與標(biāo)準(zhǔn)點(diǎn)的距離作為權(quán)函數(shù)的參數(shù)。得出權(quán)重。
3.權(quán)函數(shù)：
4.法方程的推導(dǎo)：
對于任意函數(shù) h(x) 和 g(x), 引入記號：
那么，公式4可寫成
寫成矩陣形式：寫成矩陣形式：
由上面的法方程, 解得 a(x).
然后求解得出A(x),B(x) 求解得出α(x)
5.擬合流程
6.網(wǎng)格化
這里說明一下為什么要網(wǎng)格化，網(wǎng)格化主要是選取標(biāo)準(zhǔn)點(diǎn)，并以標(biāo)準(zhǔn)點(diǎn)來劃分支撐域，確定支撐域半徑和支撐域內(nèi)的節(jié)點(diǎn) x
7.python代碼
#主題部分
X=np.arange(-0.9,0.9,0.05)
# 數(shù)據(jù)點(diǎn)x個數(shù)
M=len(x)
# 基函數(shù)個數(shù)
N=2
p=np.zeros((M,2))
Y=[]
for XX in X:
w = np.zeros((M,1))
d=0.1 # 影響區(qū)域的半徑
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)
#其他函數(shù)部分
# 權(quán)函數(shù)
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
# 權(quán)函數(shù)記號(pm,pn)的計(jì)算
def pm_pn(w,x,p,m,n):
# x為數(shù)據(jù)點(diǎn)，w為支撐域的權(quán)重，M為數(shù)據(jù)點(diǎn)個數(shù) p1,p2為傳入的數(shù)值
pmn=0
M=len(x)
# i代表數(shù)據(jù)點(diǎn),m n代表(pm,pn)的下標(biāo)
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) #數(shù)據(jù)點(diǎn)個數(shù)
m=p.shape[1] # 基函數(shù)個數(shù)
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)函數(shù)
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
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 2 11:39:26 2019
@author: yxh
"""
# Python實(shí)現(xiàn)MLS曲線擬合
import time
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#其他函數(shù)部分
# 權(quán)函數(shù)
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
# 權(quán)函數(shù)記號(pm,pn)的計(jì)算
def pm_pn(w,x,p,m,n):
# x為數(shù)據(jù)點(diǎn)，w為支撐域的權(quán)重，M為數(shù)據(jù)點(diǎn)個數(shù) p1,p2為傳入的數(shù)值
pmn=0
M=len(x)
# i代表數(shù)據(jù)點(diǎn),m n代表(pm,pn)的下標(biāo)
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) # 數(shù)據(jù)點(diǎn)個數(shù)
m=p.shape[1] # 基函數(shù)個數(shù)
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
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
path = '/home/yxh/vispy/001095_curve.txt'
pos = np.loadtxt(path)
pos_rgb = pos[:, -1]
idxs = np.where(pos_rgb == 1)
lane = pos[idxs]
x = lane[:, 0]
y = lane[:, 1]
z = lane[:, 2]
lane = lane[lane[:, 0].argsort()]
nums = 0
num = np.uint8(lane.shape[0]/50) + 1
lane_nums = np.zeros((num, 3))
for i in range(lane.shape[0]):
if(i==0 or i%50==0):
lane_nums[nums, :] = lane[i, 0:3]
nums += 1
X = lane_nums[:, 0]
#Y = lane_nums[:, 1]
Z = lane_nums[:, 2]
#X=np.arange(np.min(x, axis = 0), np.max(x, axis = 0), 1.0)
# 數(shù)據(jù)點(diǎn)x個數(shù)
M=len(x)
# 基函數(shù)個數(shù)
N=6
p=np.zeros((M,N))
Y=[]
time_begin = time.time()
ids = 0
for XX in X:
w = np.zeros((M,1))
d=4.0 # 影響區(qū)域的半徑
for i in range(0,M):
w[i]=W_fun(d, x[i], XX)
p[i][0]=1
p[i][1]=x[i]
p[i][2]=z[i]
p[i][3]=pow(x[i], 2)
p[i][4]=x[i]*z[i]
p[i][5]=pow(z[i], 2)
A=fun_A(x,w,p)
B=fun_B(y,w,p)
print len(A), len(B)
a=np.linalg.solve(A,B)
Y.append(a[0]+a[1]*XX+a[2]*Z[ids]+a[3]*pow(XX, 2)+a[4]*XX*Z[ids]+a[5]*pow(Z[ids], 2))
ids += 1
time_end = time.time()
print 'time_cost:', time_end-time_begin
Z = np.zeros((X.shape))
print Z, Z.shape
f1 = np.polyfit(X, Y, 2) # Least squares polynomial fit. use 4 can get best fitting
p1 = np.poly1d(f1)
print 'xishu:', f1
print 'fangcheng:', p1
yval = p1(X)
final_pos = np.zeros((len(Y), 3))
final_pos[:, 0] = X
final_pos[:, 1] = yval
final_pos[:, 2] = np.mean(z, axis = 0)
#Z = final_pos[:, 2]
print final_pos
print ">>>>>>>>>>>>>>>>>final_pos[0]:", final_pos[0]
print ">>>>>>>>>>>>>>>>>final_pos[-1]:", final_pos[-1]
print final_pos.shape
"""
plot1 = plt.plot(X, Y, '-s', label='original values')
plot2 = plt.plot(X, yval, 'r', label='ployfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4)
plt.title('ployfitting')
plt.show()
"""
fig = plt.Figure()
ax = plt.gca(projection='3d')
#ax3 = plt.axes(projection='3d')
ax.set_title('3D-curve')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
figure1 = ax.plot(X, Y, Z, c = 'r') #red
figure2 = ax.plot(X, yval, Z) #blue
plt.show()
#xx, yy = np.meshgrid(X, Y)
#print xx.shape
#zz = np.zeros((xx.shape))
#print zz.shape
#ax3.plot_surface(xx, yy, zz, cmap='rainbow')
#plt.show()