移动最小二乘的曲线曲面拟合

Question

LS Approximation(Least-squares)

Question: 给定 N 个点 $x_i\in\mathbb{R}^d,where\ i\in[1…N]$，希望得到 $f(x)$ 能近似出点 $x_i$ 处的 $f_i$ 值。
误差函数为 $J_{LS}=\sum_i||f(X_i)-f_i||^2$。那么得到最小化问题：

$
min_{f_x\in\prod_m^d}\sum_i||f(\mathbb{x}_i)-f_i||^2
$

$f$ 取自$\prod_m^d$,d 维空间上为 total degree 为 m 的多项式空间，则可以写作：

$$
f(x)=b(x)^Tc=b(x)\cdot c
$$

其中：$\mathbb{b}(x)=[b_1(x),…,b_k(x)]^T$是多项式的基向量，$\mathbb{c}=[c_1,…,c_k]^T$是需要最小化求解的未知系数向量。

多项式的基：

• $(m=2,d=2)，\mathbb{b(x)}=[1,x,y,x^2,xy,y^2]^T$
• $\mathbb{R^3}(m=1,d=3)$的线性拟合，$\mathbb{b(x)}=[1,x,y,z]^T$
• 任意维度常数拟合，$\mathbb{b(x)}=[1]$

Solution: 误差函数$J_{LS}$的偏导为 0.$\bigtriangledown J_{LS}=0,\bigtriangledown=[\partial/\partial_{c_1},…,\partial/\partial_{c_k}]$

因此，可以得到一个线性方程组(LSE)：

$$
\begin{align}
\partial J_{LS}/\partial c_1=0 &: \sum_i2b_1(\mathbb{x_i})[\mathbb{b(x_i)}^T\mathbb{c}-f_i]=0\\
\partial J_{LS}/\partial c_2=0 &: \sum_i2b_2(\mathbb{x_i})[\mathbb{b(x_i)}^T\mathbb{c}-f_i]=0\\
&\ \vdots\\
\partial J_{LS}/\partial c_k=0 &: \sum_i2b_k(\mathbb{x_i})[\mathbb{b(x_i)}^T\mathbb{c}-f_i]=0\\
\end{align}
$$

可以向量表示为：

$$
\begin{aligned}
& 2\sum_i\mathbb{[b(x_i)b(x_i)}^T\mathbb{c}-\mathbb{b(x_i)}f_i]=0 \\
& \Rightarrow \\
&c=[\sum_ib(x_i)b(x_i)^T]^{-1}\sum_ib(x_i)f_i
\end{aligned}
$$

Example.

在$\mathbb{R}^2$拟合一个二元二次多项式，i.e. $d=2,m=2$，那么有$\mathbb{b(x)}=[1,x,y,x^2,xy,y^2]^T$，那么线性方程组为：

$$
\begin{bmatrix}
1 & x_i & y_i & x_i^2 & x_iy_i & y_i^2 \\
x_i & x_i^2 & x_iy_i & x_i^3 & x_i^2y_i & x_iy_i^2 \\
y_i & x_iy_i & y_i^2 & x_i^2y_i & x_iy_i^2 & y_i^3 \\
x_i^2 & x_i^3 & x_i^2y_i & x_i^4 & x_i^3y_i & x_i^2y_i^2 \\
x_iy_i & x_i^2y_i & x_iy_i^2 & x_i^3y_i & x_i^2y_i^2 & x_iy_i^3 \\
y_i^2 & x_iy_i^2 & y_i^3 & x_i^2y_i^2 & x_iy_i^3 & y_i^4
\end{bmatrix}
\begin{bmatrix}
c_1\\c_2\\c_3\\c_4\\c_5\\c_6
\end{bmatrix}=
\sum_i
\begin{bmatrix}
1\\x_i\\y_i\\x_i^2\\x_iy_i\\y_i^2
\end{bmatrix}
f_i
$$

考虑一组二维平面上的点$P_i={(1,1),(1,-1),(-1,1),(-1,-1),(0,0),(1,0),(-1,0),(0,1),(0,-1)}$有两组对应的函数值$f_i^1={1.0, -0.5, 1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0}$和$f_i^2={1.0, -1.0, 0.0, 0.0, 1.0, 0.0, -1.0, -1.0, 1.0}$。

图 1 为该数据的拟合：

图 1. 二元二次多相似局部拟合：第一行为 2 组 9 个数据点，第二行是 LS 拟合函数。对应的参数向量$[c_1,…,c_6]^T$分别左为$[.0.834,.0.25,0.75,0.25,0.375,0.75]^T$，右为$[0.334,0.167,0.0,.0.5,0.5,0.0]^T$.

---

WLS Approximation(Weighted Least Squates)

在WLS方法中使用误差函数$J_{WLS}=\sum_i\theta(||\mathbb{\bar x-x_i}||)||f(\mathbb{x_i})-f_i||^2$，$\bar x \in \mathbb{R}^d$为固定点。因此需要最小化：

$$
min_{f_x\in\prod_m^d}\sum_i\theta(||\mathbb{\bar x-x_i}||)||f(\mathbb{x}_i)-f_i||^2\tag 1
$$

其中：权重 $\theta(d)$ 由数据点 $\mathbb{x_i}$ 到 $\mathbb{\bar x}$ 之间的欧氏距离表示。

未知系数也与上述权重有关，则：

$$
\begin{aligned}
f_{\mathbb{\bar x}}(\mathbb{x})&=\mathbb{b(x)^Tc(\bar x)}=\mathbb{b(x)\cdot c(\bar x)},||x-\bar x||<h \\
b(x)&=[b_1(x),…,b_k(x)]^T \\
Gaussian:\theta (d)&=e^{-\frac{d^2}{h^2}} \\
Wendland:\theta (d)&=(1-\frac{d}{h})^4(4\frac{d}{h}+1) \\
c(\bar x)&=[\sum_i\theta (d_i)b(x_i)b(x_i)^T]^{-1}\sum_i\theta (d_i)b(x_i)f_i
\end{aligned}
$$

MLS Approximation(Moving least squares)

与 WLS 很类似，不同点：

n move this point ($\bar x$) over the entire parameter domain, where a weighted least squares fit is computed and evaluated for each point individually.

因此：

$$
\begin{aligned}
f(x)&=f_x(x),min_{f_x\in\prod_m^d}\sum_i\theta(||x-x_i||)||f_x(x_i)-f_i||^2 \\
f_{\bar x}(x)&=b(x)\cdot c(\bar x) \\
c(\bar x)&=[\sum_i\theta (d_i)b(x_i)b(x_i)^T]^{-1}\sum_i\theta (d_i)b(x_i)f_i
\end{aligned}
$$

---

拟合函数的建立

在拟合区域的局部子域上，拟合函数

$$
f(x)=\sum_{i=1}^m\alpha_{i}(x)p_i(x)=p^T(x)\alpha(x)\tag{1}
$$

其中：

$\alpha(x)=[\alpha_1(x),\alpha_2(x),…,\alpha_m(x)]^T$ 为待求系数，它是坐标 x 的函数

$p(x)=[p_1(x),p_2(x),…,p_m(x)]^T$ 称为基函数，它是一个 $k$ 阶完备的多项式

$m$ 是基函数的项数

对于二维问题：

线性基：$p(x)=[1,x,y]^T，m=3$

二次基：$p(x)=[1,x,y,x^2,xy,y^2]^T，m=6$

考虑加权离散 $L_2$ 范式（向量 $\mathbf{x}=[x_1,x_2,…,x_n]$ 的 $L_2$ 范式 $L_2=(\sum_{i=1}^nx_i^2)^{1/2}$

$$
\begin{aligned}
J&=\sum_{i=1}^nw(x-x_i)[f(x)-y_i]^2 \\
&=\sum_{i=1}^nw(x-x_i)[p^T(x_i)\alpha(x)-y_i]^2
\end{aligned}\tag{3}
$$

其中：

$n$ 是影响区域内节点的数目

$f(x)$ 是拟合函数

$y_i$ 是 $x=x_i$ 处的节点值，$y_i=y(x_i)$

$w(x-x_i)$ 是节点 $x_i$ 的权函数

为确定系数 $\alpha(x)$，$J$ 应取极小值。对 $\alpha$ 求导：

$$
\frac{\partial J}{\partial \alpha}= A(x)\alpha(x)-B(x)y=0\tag{4}
$$

$$
\alpha(x)= A^{-1}(x)B(x)y\tag{5}
$$

其中：

$$
A(x)= \sum_{i=1}^nw(x-x_i)p(x_i)p^T(x_i)\tag{6} \\
$$

或者（线性二维曲线)：

$$
\begin{aligned}
A(x)&= \sum_{i=1}^np^T(x_i)w(x-x_i)p(x_i) \\
&= \begin{bmatrix}
\sum_{i=1}^nw(x-x_i) & \sum_{i=1}^nx_iw(x-x_i) \\
\sum_{i=1}^nx_iw(x-x_i) & \sum_{i=1}^nx_i^2w(x-x_i)
\end{bmatrix}
\end{aligned}
$$

$$
B(x)= [w(x-x_1)p(x_1),w(x-x_2)p(x_2),…,w(x-x_n)p(x_n)]\tag{7}
$$

$$
y^T= [y_1,y_2,…,y_n]\tag{8}
$$

把式 (5) 代入式 (1) ，就可以得到 MLS 拟合函数：

$$
f(x)=\sum_{i=1}^n\Phi_i^k(x)y_i=O^k(x)y\tag{9}
$$

其中，$O^k(x)$叫 形函数 ，$k$表示基函数的 阶数 ：

$$
O^k(x)=[\Phi_1^k,\Phi_2^k,…,\Phi_n^k]=p^T(x)A^{-1}(x)B(x)\tag{10}
$$

如果$k=0$，则基函数$p(x)={1}$，这时的形函数为 Shepard 函数：

$$
O_i^{Shepard}(x)=\frac{w(x-x_i)}{\sum_{j=1}^nw(x=x_j)}\tag{11}
$$

即使基函数$p(x)$为多项式，式 (9) 中的$f(x)$也不再是多项式。

若基函数 $p\in C^r$，权函数 $w\in C^s$，则拟合函数 $f\in C^{min(r,s)}$。

权函数

移动最小二乘中的权函数$w(x-x_i)$，其在$x$的一个子领域内不等于 0，在这个子域外全为 0，这个子域称为权函数的支持域（即$x$的影响区域)。

一般选择圆形作为权函数的支持域，其半径记为$s_{max}$

只有包含在支持域内的数据对$x$的取值有影响

权函数$w(x-x_i)$应该是非负的，并且随着$||x-x_i||_2$的增加单调递减。

权函数还应该具有一定的光滑性，因为拟合函数会继承权函数的连续性。

连续性：若权函数$w(x-x_i)$是$C^1$阶连续的，则拟合函数也是$C^1$阶连续的。

常用的权函数是样条函数，记$s=x-x_i，\bar s=\frac{s}{s_{max}}$则三次样条函数如式 (12) 所示：

$$
\begin{aligned}
w (\bar s) &= \frac{2}{3}-4\bar s^2 + 4\bar s^3 & (\bar s\le\frac{1}{2}) \\
w (\bar s) &= \frac{4}{3}-4\bar s+4\bar s^2-\frac{4}{3}\bar s^3 & (\frac{1}{2}<\bar s \le 1) \\
w (\bar s) &= 0 & (\bar s>1)
\end{aligned}
\tag{12}
$$

支持域应该包含足够多的节点，使得式 (6) 中$A(x)$可逆。

如果式 (1) 使用的是线性基函数，则 曲线拟合 的支持域内应该至少包含不重叠的 2 个节点，曲面拟合 的支持域内应该至少包含不在同一条直线上的 3 个节点

MLS 拟合流程

基本思想：先讲拟合区域网格化，然后用公式 (9) 求出网格上节点值，最后连接网格节点形成拟合曲线（曲面）。

• 将拟合区域网格化；
• 对每个网格点 x 进行循环：
  • 确定网格点 x 的支持域（半径) 的大小；
  • 确定包含在 x 的影响区域内的节点；
  • 计算行函数$O^k(x)$；
  • 计算网格点 x 处的函数值
• 结束网格点循环；
• 连接网格点形成拟合曲线（曲面)。

误差分析

为分析曲面拟合误差，定义误差项：

$$
error=\int_U[f^{EXACT}(x,y)-f^{MLS}(x,y)]^2dU\tag{14}
$$

其中，$f^{EXACT}(x,y)$是精确值，$f^{MLS}(x,y)$是移动最小二乘的计算值。

---

Applications

The MLS surface of a point-set with varying density (the
density is reduced along the vertical axis from top to bottom). The
surface is obtained by applying the projection operation described
by Alexa et. al. [2003]. Image courtesy of Marc Alexa.

---

Coding

• Matlab 原理版本
• C++版本