样条曲线长度--数值积分

转载自：https://www.cnblogs.com/flytrace/p/8413255.html

文章目录

1.根据长度求坐标思路
2.用数值积分实现s的积分
3.牛顿迭代求解s对应的t
4.直接推导二阶贝塞尔曲线的长度积分解析式

前几天在做摄像机轨道时，解决匀速定长运动问题，刚好理解了题述问题。

首先做完这个曲线上匀速运动的编程实现后，是十分愉悦的。这个问题里，两个科学史上的伟人，牛顿和高斯，相继大显身手，怎能不让人膜拜？我之前并不了解所涉及的数值分析方法，所以觉得很神奇。

以下为简洁我会使用向量的表示方法，使用大写数字代表向量。以catmullrom样条曲线(spline)为例，对于一个关于t的多次样条曲线:

Q(t)=At3+Bt2+Ct+DQ(t) = At^{3} + Bt^{2} + Ct + D Q(t)=At3+Bt2+Ct+D

1.根据长度求坐标思路

一般来说实际中，ttt的范围都是[0, 1]。当ttt在[0, 1]中滑动时，Q(t)Q(t)Q(t)的集合即为一条平滑曲线。注意到当ttt以匀速移动时，Q(t)Q(t)Q(t)在沿曲线方向上并不是匀速移动。要做到这一点，显然曲线长度sss是需要考虑的。假设在ttt点处我们有曲线长度，（根据sss长度，计算控制量ttt，再计算坐标Q(t)Q(t)Q(t)）。
s=L(t)s = L(t) s=L(t)
那么s长度所对应的t即为上面方程的反函数:
t=L−1(s)t = L^{-1}(s) t=L−1(s)
反带入Q(t)Q(t)Q(t)，得
Q(t)=Q(L−1(s))=Z(s)Q(t) = Q(L^{-1}(s)) = Z(s) Q(t)=Q(L−1(s))=Z(s)

sss的区间显然是[0, 曲线总长度]，于是我们可以对sss做匀速递增，从而得到Q(t)Q(t)Q(t)的位置，从而完成匀速插值。

难点在于求S=L(t)S=L(t)S=L(t)，这实际上是一个求积分的问题。很遗憾，一般来说很难得到一个解析的积分公式的，实际中是通过数值计算方法得到积分值。

sss通过数值方法得到，也很难求逆得到ttt，可以通过牛顿迭代，或二分查找得到。

2.用数值积分实现s的积分

数值积分使用Guass-lengendre积分来实现，通过选定的几个积分点及权重，它神奇而快速的得到精度很高的积分值。我们直接来重点吧，该方法需要一张n阶权重表，比如5阶权重表:

及一个计算公式:
∫abf(x)dx=b−a2∑i=1nwif(b−a2xi+b+a2)\int_{a}^{b}f(x)dx = \frac{b-a}{2}\sum_{i=1}^{n}{w_{i}f(\frac{b-a}{2}x_{i}+\frac{b+a}{2})} ∫abf(x)dx=2b−ai=1∑nwif(2b−axi+2b+a)
对于样条曲线，我们有:
dsdt=∣∣Q′(t)∣∣ds=∣∣Q′(t)∣∣dts=∫0t∣∣Q′(x)∣∣dx\frac{ds}{dt} = ||Q'(t)|| \\ ds = ||Q'(t)||dt \\ s = \int_{0}^{t}||Q'(x)||dx dtds=∣∣Q′(t)∣∣ds=∣∣Q′(t)∣∣dts=∫0t∣∣Q′(x)∣∣dx

Q′(t)Q'(t)Q′(t)为Q(t)Q(t)Q(t)的导
Q′(t)=3At2+2Bt+CQ'(t) = 3At^{2} + 2Bt + C Q′(t)=3At2+2Bt+C
∣∣Q′(t)∣∣||Q'(t)||∣∣Q′(t)∣∣为该点处的长度:
∣∣Q′(t)∣∣=Qx′(t)2+Qy′(t)2||Q'(t)|| = \sqrt{Q_x'(t)^2 + Q_y'(t)^2} ∣∣Q′(t)∣∣=Qx′(t)2+Qy′(t)2
具体对于我们的公式，有:
s=L(t)=∫0t∣∣Q′(x)∣∣dx=t2∑i=1nwi∣∣Q′(t2xi+t2)∣∣s= L(t) = \int_{0}^{t}||Q'(x)||dx = \frac{t}{2}\sum_{i=1}^{n}{w_{i}||Q'(\frac{t}{2}x_{i}+\frac{t}{2})||} s=L(t)=∫0t∣∣Q′(x)∣∣dx=2ti=1∑nwi∣∣Q′(2txi+2t)∣∣

将表中的值依次带入，迭代计算即可得到ttt点处的积分值sss。比我之前快速实现时所采用的线性分割方法快了100倍有没有？特别是该方程的普适性，对于多阶多项式都适用而且精度非常高。有种魔性在这里，高斯是神。

3.牛顿迭代求解s对应的t

现在s=L(t)s=L(t)s=L(t)对我们来说是已知了，那么如何求sss所对应的ttt呢?这实际上是求方程:
L(t)−s=0L(t) - s = 0 L(t)−s=0
的根。因为L(t)L(t)L(t)没有解析表达，我们将使用牛顿迭代法:
b=a−f(a)f′(a)b= a - \frac{f(a)}{f'(a)} b=a−f′(a)f(a)
将得到的b的值作为a再带入上述公式，迭代几次即可逼近s=f(x)s = f(x)s=f(x)的根。具体到我们的公式：
b=a−L(a)−sL′(a)b= a - \frac{L(a) - s}{L'(a)} b=a−L′(a)L(a)−s
对于特定的sss，我们首先要拟定一个初值aaa，因为样条曲线性质良好，而且一般说来实际中不会允许2点间出现陡峭的曲线，一定会多加一些点令2点间曲线形状为凸或凹的，所以我们可以近似认为ttt在区间[0, 1]的比例，近似等于sss与总长度的比例。而显然总长为L(1),通过前面的高斯积分已经可以算出，于是初始值aaa我们可拟定为：
a=sL(1)a = \frac{s}{L(1)} a=L(1)s
反复迭代(3-4次的精度就不错了)，得出最终ttt值，再把该值反带入Q(t)Q(t)Q(t)中，得到对于给定s值所对应的位置，至此我们完成了样条曲线的重参数化。

4.直接推导二阶贝塞尔曲线的长度积分解析式

https://blog.csdn.net/LANGZI7758521/article/details/52101672