移動最小二乘法在點雲平滑和重採樣中的應用

時間 2019-12-10

原文原文鏈接

目錄php

導言

這篇博客的目的是介紹移動最小二乘法（moving least squares）在點雲平滑和重採樣的應用。開頭先簡要介紹最小二乘法，用兩種不一樣的方法來求解最小二乘法，並給出一個具體的算例、代碼、數據。目前關於最小二乘法的博客和網上的討論有不少，我不必重複作這個工做（實際上是我太菜不能形象的講出來哈），我想作的是提供一個簡要的歸納，以及給出一個具體的例子幫助讀者理解。接下來，介紹加權最小二乘法，在一樣的數據上面，跑加權最小二乘法，看看擬合的結果又如何。以後，簡略分析一下PCL是如何實現的點雲平滑和重採樣。參考的文獻列在最後，比較有價值的是[1][2]，基本上是看這兩個理解的移動最小二乘法。因爲我學力尚淺，自知這篇文章有不少的不足、錯誤，但願讀者不吝指出。歡迎理性的討論和交流。html

最小二乘法

對於方程組\(A\mathbf{x}=b\)，咱們要去尋找一個\(Ax'\)最佳逼近b，即有
\[|b - Ax'| \le |b - Ax|\]
抽象點說，在\(A\)的列空間中找到最逼近b的那個點。python

使用投影的方法來求解

將向量\(b\)投影到\(A\)的列空間獲得\(b'\)，\((b - b')\)是正交於A的列空間，即有
\[A^T*(b - Ax') = 0\]
若是\((A^T*A)\)可逆：
\[x' = (A^T*A)^{-1}*A^T*b\]數組

例1

這組數據，來自於[2]。畫出散點圖。dom

x = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
y = [0, 4, 5, 14, 15, 14.5, 14, 12, 10, 5, 4]

首先考慮使用\(y = mx + n\)來擬合（這裏用m,n而不用a,b是怕讀者混淆）。須要解決的問題是數據如何裝填。A矩陣放什麼，b放什麼，x又是什麼。
欲求的是\(m,n\)。\(x = [m \quad n]^T\)，\(b\)放入的是每個因變量。欲求\(Ax = b\)，那麼每一行A乘以x，會等於對應的每一行b:
\[A = \begin{bmatrix}x_1&1\\x_2&1\\...&...\\x_n&1\\\end{bmatrix}\]
\[b = \begin{bmatrix}y_1\\y_2\\...\\y_n\end{bmatrix}\]函數

A = np.empty((0, 2))
B = np.empty((0, 1))
for (a, b) in zip(x, y):
    row1 = np.array([a, 1])
    row2 = np.array([b])
    A = np.vstack([A, row1])
    B = np.vstack([B, row2])
ans = np.linalg.inv(np.mat(A.transpose()) * np.mat(A)) * np.mat(A.transpose()) * np.mat(B)
print(ans)
xx = np.arange(0, 1, 0.01)
yy = ans[0, 0] * xx + ans[1, 0]
plt.plot(xx, yy)
plt.scatter(x, y, c='r')
plt.show()

擬合結果：
spa

例2

但是這個擬合的結果太差勁了吧！咱們觀察數據，以爲這個能夠用二次函數去擬合，那麼可使用新的擬合函數：\(y = ax^2 + bx + c\)3d

\[A = \begin{bmatrix}x_1^2&x_1&1\\x_2^2&x_2&1\\...&...\\x_n^2&x_n&1\\\end{bmatrix}\]
\[x = [a\quad b\quad c]^T\]code

再擬合一次：此次擬合的還行。\(a = -53.96270396, b = 57.05361305, c = -0.77622378\)orm

A = np.empty((0, 3))
B = np.empty((0, 1))
for (a, b) in zip(x, y):
    row1 = np.array([a*a, a, 1])
    row2 = np.array([b])
    A = np.vstack([A, row1])
    B = np.vstack([B, row2])
print(A)
print(B)
ans = np.linalg.inv(np.mat(A.transpose()) * np.mat(A)) * np.mat(A.transpose()) * np.mat(B)
print(ans)
xx = np.arange(0, 1, 0.01)
yy = ans[0, 0] * xx**2 + ans[1, 0] * xx + ans[2, 0]
plt.plot(xx, yy)
plt.scatter(x, y, c='r')
plt.show()

使用求偏導的方法來求解

另外一種理解最小二乘法的方法是，找到使因變量的偏差平方和最小的參數來求解。假設咱們有二維數據點集\((x_i, y_i)\)，對他擬合的函數爲\(f(x)\)。
\[J_{LS} = \sum^{k}_{n=1}(f(x_i) - y_i)^2\]
因而問題化解爲求\(min(J_{LS})\)。
爲何使用偏差平方和呢？我以爲溝通這個方法和上一個方法的關鍵在於距離！偏差平方和，實際上是因變量的向量\((y_1, y_2, ... , y_n)\)和擬合函數對應函數值的向量\((f(x_1), f(x_2), ... , f(x_n))\)的距離平方。要使到這個距離最小，那就是投影！其實這兩種方法本質上都是利用投影來使到誤差最小。不一樣之處在於一種使用了內積爲0的方式來求，一種使用求偏導算最小值的方式來求。

[1][2]都是用了一個係數向量\(\alpha\)，和基函數\(p(x)\)來表示\(f(x) = p(x)^T * \alpha\)。

若是有一個自變量x，最高階爲3，那麼\(p(x) = [1,x,x^2,x^3]^T\)
若是有兩個自變量x、y，最高階爲2，那麼\(p(x) = [1, x, y, x^2, xy, y^2]^T\)

係數向量是咱們想要求的向量, \(\alpha = [c_1, c_2, ..., c_n]\)。

對\(J_{LS}\)分別求\(c_1, c_2, ..., c_n\)的偏導，具體過程參考[1]
\[\frac{\partial J_{LS}}{\partial \alpha} = \sum^n_{i=1} 2*p(x_i)*(p(x_i)^T * \alpha - y_i) \]

令偏導爲0
\[\alpha = [\sum^n_{i=1}(p(x_i) * p(x_i)^T)]^{-1} * \sum^n_{i=1}p(x_i)*y_i\]

注意到\(p(x_i)\)是一個\(n\times 1\)的向量，他的轉置則是\(1\times n\)，因而它和它的轉置相乘是一個\(n\times n\)的矩陣。

例子3

使用例子2的函數去擬合：\(y = a + bx + cx^2\)
基函數\(p(x) = [1, x, x^2]^T\), 係數向量\(\alpha = [a, b, c]\)
\[p(x) * p(x)^T = \begin{bmatrix}1&x&x^2\\x&x^2&x^3\\x^2&x^3&x^4\\\end{bmatrix}\]
擬合出來的結果爲：\(a=-0.77622378,b=57.05361305,c=-53.96270396\)。異曲同工啊！

sumx = sumxx = sumxxx = sumxxxx = sumf = sumxf = sumxxf = 0
for (a, b) in zip(x, y):
    sumx += a
    sumxx += a ** 2
    sumxxx += a ** 3
    sumxxxx += a ** 4
    sumf += b
    sumxf += a * b
    sumxxf += a * a * b
A = np.array([[len(x), sumx, sumxx],
              [sumx, sumxx, sumxxx],
              [sumxx, sumxxx, sumxxxx]])
B = np.array([sumf, sumxf, sumxxf])
ans = np.linalg.solve(A, B)
print(ans)
xx = np.arange(0, 1, 0.01)
yy = ans[0] + ans[1] * xx + ans[2] * xx**2
plt.plot(xx, yy)
plt.scatter(x, y, c='r')
plt.show()

加權最小二乘法

[2]論文中：

引入緊支（ Compact Support）概念，認爲點 x 處的值 y 只受 x 附近子域內節點影響

因變量只受到自變量某一鄰域影響，引入一個權重函數，給重要的地方加權，不重要的地方削弱它對因變量的影響。

\[J_{WLS} = \sum^n_{i=1}\theta(s) (f(x_i) - y_i)^2\]
其中\(s\)爲想要求的自變量\(x\)到附近樣本自變量的歐幾里得距離。

若是咱們想根據一組數據去估計某一個自變量\(x\)的函數值，咱們就得計算一次！這種方法計算量相對而言大不少。

一樣這是一個最小化問題\(min(J_{WLS})\)，對它求偏導，令偏導爲0。

\[\alpha = [\sum^n_{i=1}\theta(s)(p(x_i) * p(x_i)^T)]^{-1} * \sum^n_{i=1}\theta(s)p(x_i)*y_i\]

例子4

仍然是那一組數據。使用[2]建議的加權函數來求解。只不過此次使用的是一個一次函數來擬合局部。\(p(x) = [1,x]^T\)
所謂的移動最小二乘法，找到一個函數，這個函數有一系列的局部的函數組成\(f(x) \in f_X(x)\)，本質上就是WLS。

def w(dis):
    dis = dis / 0.3
    if dis < 0:
        return 0
    elif dis <= 0.5:
        return 2/3 - 4 * dis**2 + 4 * dis**3
    elif dis <= 1:
        return 4/3 - 4 * dis + 4 * dis**2 - 4/3 * dis**3
    else:
        return 0


def mls(x_):
    sumxx = sumx = sumxf = sumf = sumw = 0
    for (a, b) in zip(x, y):
        weight = w(abs(x_ - a))
        sumw += weight
        sumx += a * weight
        sumxx += a * a * weight
        sumf += b * weight
        sumxf += a * b * weight
    A = np.array([[sumw, sumx],
                  [sumx, sumxx]])
    B = np.array([sumf, sumxf])
    ans = np.linalg.solve(A, B)
    return ans[0] + ans[1] * x_


x = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
y = [0, 4, 5, 14, 15, 14.5, 14, 12, 10, 5, 4]
xx = np.arange(0, 1, 0.01)
yy = [mls(xi) for xi in xx]
plt.plot(xx, yy)
plt.scatter(x, y, c='r')
plt.show()