转载自:https://www.cnblogs.com/flytrace/p/8413255.html
前几天在做摄像机轨道时,解决匀速定长运动问题,刚好理解了题述问题。
首先做完这个曲线上匀速运动的编程实现后,是十分愉悦的。这个问题里,两个科学史上的伟人,牛顿和高斯,相继大显身手,怎能不让人膜拜?我之前并不了解所涉及的数值分析方法,所以觉得很神奇。
以下为简洁我会使用向量的表示方法,使用大写数字代表向量。以catmullrom样条曲线(spline)为例,对于一个关于t的多次样条曲线:
Q ( t ) = A t 3 + B t 2 + C t + D Q(t) = At^{3} + Bt^{2} + Ct + D Q(t)=At3+Bt2+Ct+D
1.根据长度求坐标思路
一般来说实际中,
t
t
t的范围都是[0, 1]。当
t
t
t在[0, 1]中滑动时,
Q
(
t
)
Q(t)
Q(t)的集合即为一条平滑曲线。注意到当
t
t
t以匀速移动时,
Q
(
t
)
Q(t)
Q(t)在沿曲线方向上并不是匀速移动。要做到这一点,显然曲线长度
s
s
s是需要考虑的。假设在
t
t
t点处我们有曲线长度,(根据
s
s
s长度,计算控制量
t
t
t,再计算坐标
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)
s s s的区间显然是[0, 曲线总长度],于是我们可以对 s s s做匀速递增,从而得到 Q ( t ) Q(t) Q(t)的位置, 从而完成匀速插值。
难点在于求 S = L ( t ) S=L(t) S=L(t),这实际上是一个求积分的问题。很遗憾,一般来说很难得到一个解析的积分公式的,实际中是通过数值计算方法得到积分值。
s s s通过数值方法得到,也很难求逆得到 t t t,可以通过牛顿迭代,或二分查找得到。
2.用数值积分实现s的积分
数值积分使用Guass-lengendre积分来实现,通过选定的几个积分点及权重,它神奇而快速的得到精度很高的积分值。我们直接来重点吧,该方法需要一张n阶权重表,比如5阶权重表:
及一个计算公式:
∫
a
b
f
(
x
)
d
x
=
b
−
a
2
∑
i
=
1
n
w
i
f
(
b
−
a
2
x
i
+
b
+
a
2
)
\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)
对于样条曲线,我们有:
d
s
d
t
=
∣
∣
Q
′
(
t
)
∣
∣
d
s
=
∣
∣
Q
′
(
t
)
∣
∣
d
t
s
=
∫
0
t
∣
∣
Q
′
(
x
)
∣
∣
d
x
\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 ) = 3 A t 2 + 2 B t + C Q'(t) = 3At^{2} + 2Bt + C Q′(t)=3At2+2Bt+C -
∣
∣
Q
′
(
t
)
∣
∣
||Q'(t)||
∣∣Q′(t)∣∣为该点处的长度:
∣ ∣ Q ′ ( t ) ∣ ∣ = Q x ′ ( t ) 2 + Q y ′ ( t ) 2 ||Q'(t)|| = \sqrt{Q_x'(t)^2 + Q_y'(t)^2} ∣∣Q′(t)∣∣=Qx′(t)2+Qy′(t)2 - 具体对于我们的公式,有:
s = L ( t ) = ∫ 0 t ∣ ∣ Q ′ ( x ) ∣ ∣ d x = t 2 ∑ i = 1 n w i ∣ ∣ Q ′ ( t 2 x i + t 2 ) ∣ ∣ 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)∣∣
将表中的值依次带入,迭代计算即可得到 t t t点处的积分值 s s s。比我之前快速实现时所采用的线性分割方法快了100倍有没有?特别是该方程的普适性,对于多阶多项式都适用而且精度非常高。有种魔性在这里,高斯是神。
3.牛顿迭代求解s对应的t
现在
s
=
L
(
t
)
s=L(t)
s=L(t)对我们来说是已知了,那么如何求
s
s
s所对应的
t
t
t呢?这实际上是求方程:
L
(
t
)
−
s
=
0
L(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
)
−
s
L
′
(
a
)
b= a - \frac{L(a) - s}{L'(a)}
b=a−L′(a)L(a)−s
对于特定的
s
s
s,我们首先要拟定一个初值
a
a
a,因为样条曲线性质良好,而且一般说来实际中不会允许2点间出现陡峭的曲线,一定会多加一些点令2点间曲线形状为凸或凹的,所以我们可以近似认为
t
t
t在区间[0, 1]的比例,近似等于
s
s
s与总长度的比例。而显然总长为L(1),通过前面的高斯积分已经可以算出,于是初始值
a
a
a我们可拟定为:
a
=
s
L
(
1
)
a = \frac{s}{L(1)}
a=L(1)s
反复迭代(3-4次的精度就不错了),得出最终
t
t
t值,再把该值反带入
Q
(
t
)
Q(t)
Q(t)中,得到对于给定s值所对应的位置,至此我们完成了样条曲线的重参数化。
4.直接推导二阶贝塞尔曲线的长度积分解析式
https://blog.youkuaiyun.com/LANGZI7758521/article/details/52101672