【CAD/CG笔记】求控制点对应在B-Spline上的参数估计
求控制点Pi对应在B-Spline上的参数估计
由得到的参数估计值作为迭代初始值进行牛顿迭代（泰勒展开或者分式近似），实验能否收敛
示例：Bezier曲线控制点的参数估计
首先对于Bezier曲线上的控制点，其参数估计为
$$
P_i = \frac {i}{n}
$$
注意到求和式有下述恒等关系
$$
\sum_{i=0}^{n} \frac{i}{n} B_i(t) = \sum_{i=0}^{n} \frac{i}{n} \binom{n}{i} t^i (1 - t)^{n - i} = t
$$
证明如下：
$$
\sum_{i=0}^{n} \frac{i}{n} B_i(t) = \sum_{i=0}^{n} \frac{i}{n} \binom{n}{i} t^i (1 - t)^{n - i}
$$
其中 $B_i(t)$ 是 Bernstein 基函数：
$$
B_i(t) = \binom{n}{i} t^i (1-t)^{n-i}
$$
注意到当 ( i = 0 ) 时，该项为 0，因此求和可从 ( i = 1 ) 开始：
$$
= \sum_{i=1}^{n} \frac{i}{n} \binom{n}{i} t^i (1 - t)^{n - i}
$$
利用组合恒等式：
$$
\frac{i}{n} \binom{n}{i} = \binom{n-1}{i-1}
$$
代入后得到：
$$
= \sum_{i=1}^{n} \binom{n-1}{i-1} t^i (1 - t)^{n - i}
$$
令 j = i - 1，则求和变为：
$$
= \sum_{j=0}^{n-1} \binom{n-1}{j} t^{j+1} (1 - t)^{(n-1) - j} = t \sum_{j=0}^{n-1} \binom{n-1}{j} t^{j} (1 - t)^{(n-1) - j}
$$
根据二项式定理，求和部分为$(t + (1 - t))^{n-1} = 1^{n-1} = 1$，因此：
$$
= t \cdot 1 = t
$$
最终结果为：
$$
\sum_{i=0}^{n} \frac{i}{n} B_i(t) = \sum_{i=0}^{n} \frac{i}{n} \binom{n}{i} t^i (1 - t)^{n - i} = t
$$
附.关于组合恒等式
从左边出发：
$$
\frac{i}{n}\binom{n}{i} = \frac{i}{n} \cdot \frac{n!}{i!(n-i)!}
$$
将分子中的 $n!$ 展开为 $n \cdot (n-1)!$：
$$
= \frac{i}{n} \cdot \frac{n \cdot (n-1)!}{i!(n-i)!}
$$
约分 $n$ 和 $i$（同理 $i! = i \cdot (i-1)!$）：
$$
= \frac{(n-1)!}{(i-1)!(n-i)!}
$$
这正是组合数 $\binom{n-1}{i-1}$ 的定义：
$$
= \binom{n-1}{i-1}
$$
同理可以求出B-Spline的参数估计
假设$P_i$的参数是$r_i$，对于下面的式子
$$
\sum_{i=0}^{n} r_i B_{i,j}(t)
$$
令其尽可能的等于$t$，即
$$
\sum_{i=0}^{n} r_i B_{i,j}(t) = t
$$
那么这样得到的$r_i$就是控制点$P_i$的参数估计值
其中，B样条基函数由考克斯-德布尔递归公式定义，
$$
B_{i,0}(t) =
\begin{cases}
1, & t_i \leq t < t_{i+1} \
0, & \text{otherwise}
\end{cases}
$$
$$
B_{i,j}(t) = \frac{(t - t_i) B_{i,j-1}(t)}{t_{i+j-1} - t_i} + \frac{(t_{i+j} - t) B_{i+1,j-1}(t)}{t_{i+j} - t_{i+1}}, \quad 1 \leq j \leq n
$$
其中，约定 $0/0 = 0$；控制点的个数为 $n+1$；节点个数为 $n+j+2$；基函数的阶次为 $j+1$，次数为 $j$。
B-Spline的数学表达式为：
$$
X(t) = \sum_{i=0}^{n} B_{i,j}(t) P_i, \quad t \in [t_0, t_1], \quad 1 \leq j \leq n
$$
求B-Spline参数估计的具体做法
通过最小二乘法，我们希望使下述积分最小：
$$
\int_{0}^{1}(\sum_{i=0}^{n} r_i B_{i,j}(t) - t)^2 dx
$$
求出上述积分，展开平方项：
$$
I = \int_0^1 \left( \sum_{i=0}^{n} r_i B_{i,j}(t) - t \right)^2 dt
$$
$$
= \int_0^1 \left( \sum_{i=0}^{n} r_i B_{i,j}(t) \right)^2 dt - 2 \int_0^1 t \sum_{i=0}^{n} r_i B_{i,j}(t) dt + \int_0^1 t^2 dt
$$
记：
$$
M_{i,k} = \int_0^1 B_{i,j}(t) B_{k,j}(t) dt, \quad v_i = \int_0^1 t B_{i,j}(t) dt
$$
那么积分变为矩阵形式：
$$
I = \sum_{i,k} r_i r_k M_{i,k} - 2 \sum_i r_i v_i + \frac{1}{3}
$$
为了最小化 $I$，对 $r_i$ 求偏导：
$$
\frac{\partial I}{\partial r_m} = 2 \sum_k r_k M_{m,k} - 2v_m
$$
令其等于 0：
$$
\sum_k M_{m,k} r_k = v_m
$$
这个方程组可以表示为：
$$
Mr = v
$$
其中：
- $M$ 是 $(n+1) \times (n+1)$ 的矩阵，存储了基函数的内积，对称且半正定
- $r$ 是我们要解的参数估计向量
- $v$ 是右端向量，存储了基函数与 $t$ 的积分
最终，我们只需要求解线性方程组：
$$
r = M^{-1}v
$$
这个方程组可以用 矩阵求逆 或 LU 分解 来求解。在 MATLAB 或 Python 中，通常使用 $M \ v$ 来直接求解：
$$
r = M \backslash v
$$这是最小二乘意义下最优的参数估计值 $r_i$。
附.（22）式的详细推导
$$
I = \int_0^1 \left( \sum_{i=0}^{n} r_i B_{i,j}(t) \right)^2 dt
$$
由于平方可以写成自身与自身的乘积：
$$
I = \int_0^1 \left( \sum_{i=0}^{n} r_i B_{i,j}(t) \right) \cdot \left( \sum_{k=0}^{n} r_k B_{k,j}(t) \right) dt
$$
现在我们展开括号：
$$
I = \int_0^1 \sum_{i=0}^{n} \sum_{k=0}^{n} r_i r_k B_{i,j}(t) B_{k,j}(t) dt
$$
由于积分运算对求和符号是线性的，我们可以交换求和与积分的顺序：
$$
I = \sum_{i=0}^{n} \sum_{k=0}^{n} r_i r_k \int_0^1 B_{i,j}(t) B_{k,j}(t) dt
$$
这里，积分部分是一个 B 样条基函数的内积，我们定义它为矩阵 $M$ 的元素：
$$
M_{i,k} = \int_0^1 B_{i,j}(t) B_{k,j}(t) dt
$$
第一项就可以写成矩阵形式：
$$
\sum_{i=0}^{n} \sum_{k=0}^{n} r_i r_k M_{i,k} = \mathbf{r}^T M \mathbf{r}
$$
Matlab代码
function r = analytic_solution(j, knots) |
C++代码
|
结论
设三次 B 样条曲线，节点向量为$[0,0,0,0,u_1,u_2,…,u_n,1,1,1,1]$，
- 则控制点对应关系为： $[0,\frac{1}{3}u_1, u_1,u_2,…..,u_n,1-\frac{1-u_n}{3},1]$
设四次 B 样条曲线，节点向量为$[0,0,0,0,0,u_1,u_2,…,u_n,1,1,1,1,1]$，
- 则控制点对应关系为： $[0,\frac{u_1+u_2}{12}, \frac{u_1+u_2}{4}, \frac{u_1+u_2}{2},\frac{u_2+u_3}{2},…..,\frac{u_{n-1}+u_n}{2},1-\frac{u_1+u_2}{4} ,1-\frac{u_1+u_2}{12} ,1]$
