四元数的插值

转载：http://blog.****.net/jin_syuct/article/details/49785541

今天我们和大家分享的是四元数的插值。这里的插值指的是球面线性插值。例如，我们要模拟一下地球绕着太阳，从P1到P2。这中间的每一个位置都要用球面线性插值来做。

首先我们聊一聊线性插值：

x=x1-x2,t是插值系数，则lerp(x1,x2,t)=x1+t*x表示x1到x2的插值。

四元数的插值

q=(p-1)*p1,插值系数为t，则p到p1的插值为：slerp(p,p1,t)=p*((p-1)p1)t，表示p的逆乘以p1，他们乘积的t次方，乘以p

但是上面的公式在编程的时候用的挺麻烦的，所以我们使用下面的公式：

旋转插值，想开头所说的，这里说的插值是球面插值，是在3D空间中旋转的，因此我们可以将它等价于旋转插值。

例如：向量V1，V0

W=v1-v0

vt=v0+tw

通过上图可知，VT=K0V0+K1V1.在这里，V1，V0，VT都是单位向量，V1和K1V1平行。

因此，我们可以求得sinw=sintw/k1,则k1=sintw/sinw,同样可得 k0=sin(1-t)w/sinw.

所以代入VT，就可以求出VT了。

四元数也是一个道理，注意slerp这个函数返回的是一个四元数。

以上就是四元数插值的公式，在编程实现的时候我们要注意

使用点成来纠结夹角W，还有就是放夹角W很小的时候，sinw很小，但是cosw趋于1，因此，就变成了线性插值。

在这里，我们贴出cocos2dx 3.6中四元数的实现代码，大家学习学习。

[cpp] view plain copy

print?

1. /** 
2.  * Interpolates between two quaternions using spherical linear interpolation. 
3.  * 
4.  * Spherical linear interpolation provides smooth transitions between different 
5.  * orientations and is often useful for animating models or cameras in 3D. 
6.  * 
7.  * Note: For accurate interpolation, the input quaternions must be at (or close to) unit length. 
8.  * This method does not automatically normalize the input quaternions, so it is up to the 
9.  * caller to ensure they call normalize beforehand, if necessary. 
10.  * 
11.  * @param q1x The x component of the first quaternion. 
12.  * @param q1y The y component of the first quaternion. 
13.  * @param q1z The z component of the first quaternion. 
14.  * @param q1w The w component of the first quaternion. 
15.  * @param q2x The x component of the second quaternion. 
16.  * @param q2y The y component of the second quaternion. 
17.  * @param q2z The z component of the second quaternion. 
18.  * @param q2w The w component of the second quaternion. 
19.  * @param t The interpolation coefficient. 
20.  * @param dstx A pointer to store the x component of the slerp in. 
21.  * @param dsty A pointer to store the y component of the slerp in. 
22.  * @param dstz A pointer to store the z component of the slerp in. 
23.  * @param dstw A pointer to store the w component of the slerp in. 
24.  */  
25. static void slerp(float q1x, float q1y, float q1z, float q1w, float q2x, float q2y, float q2z, float q2w, float t, float* dstx, float* dsty, float* dstz, float* dstw);  

[cpp] view plain copy

print?

1. void Quaternion::slerp(float q1x, float q1y, float q1z, float q1w, float q2x, float q2y, float q2z, float q2w, float t, float* dstx, float* dsty, float* dstz, float* dstw)  
2. {  
3.     // Fast slerp implementation by kwhatmough:  
4.     // It contains no division operations, no trig, no inverse trig  
5.     // and no sqrt. Not only does this code tolerate small constraint  
6.     // errors in the input quaternions, it actually corrects for them.  
7.     GP_ASSERT(dstx && dsty && dstz && dstw);  
8.     GP_ASSERT(!(t < 0.0f || t > 1.0f));  
9.   
10.     if (t == 0.0f)  
11.     {  
12.         *dstx = q1x;  
13.         *dsty = q1y;  
14.         *dstz = q1z;  
15.         *dstw = q1w;  
16.         return;  
17.     }  
18.     else if (t == 1.0f)  
19.     {  
20.         *dstx = q2x;  
21.         *dsty = q2y;  
22.         *dstz = q2z;  
23.         *dstw = q2w;  
24.         return;  
25.     }  
26.   
27.     if (q1x == q2x && q1y == q2y && q1z == q2z && q1w == q2w)  
28.     {  
29.         *dstx = q1x;  
30.         *dsty = q1y;  
31.         *dstz = q1z;  
32.         *dstw = q1w;  
33.         return;  
34.     }  
35.   
36.     float halfY, alpha, beta;  
37.     float u, f1, f2a, f2b;  
38.     float ratio1, ratio2;  
39.     float halfSecHalfTheta, versHalfTheta;  
40.     float sqNotU, sqU;  
41.   
42.     float cosTheta = q1w * q2w + q1x * q2x + q1y * q2y + q1z * q2z;  
43.   
44.     // As usual in all slerp implementations, we fold theta.  
45.     alpha = cosTheta >= 0 ? 1.0f : -1.0f;  
46.     halfY = 1.0f + alpha * cosTheta;  
47.   
48.     // Here we bisect the interval, so we need to fold t as well.  
49.     f2b = t - 0.5f;  
50.     u = f2b >= 0 ? f2b : -f2b;  
51.     f2a = u - f2b;  
52.     f2b += u;  
53.     u += u;  
54.     f1 = 1.0f - u;  
55.   
56.     // One iteration of Newton to get 1-cos(theta / 2) to good accuracy.  
57.     halfSecHalfTheta = 1.09f - (0.476537f - 0.0903321f * halfY) * halfY;  
58.     halfSecHalfTheta *= 1.5f - halfY * halfSecHalfTheta * halfSecHalfTheta;  
59.     versHalfTheta = 1.0f - halfY * halfSecHalfTheta;  
60.   
61.     // Evaluate series expansions of the coefficients.  
62.     sqNotU = f1 * f1;  
63.     ratio2 = 0.0000440917108f * versHalfTheta;  
64.     ratio1 = -0.00158730159f + (sqNotU - 16.0f) * ratio2;  
65.     ratio1 = 0.0333333333f + ratio1 * (sqNotU - 9.0f) * versHalfTheta;  
66.     ratio1 = -0.333333333f + ratio1 * (sqNotU - 4.0f) * versHalfTheta;  
67.     ratio1 = 1.0f + ratio1 * (sqNotU - 1.0f) * versHalfTheta;  
68.   
69.     sqU = u * u;  
70.     ratio2 = -0.00158730159f + (sqU - 16.0f) * ratio2;  
71.     ratio2 = 0.0333333333f + ratio2 * (sqU - 9.0f) * versHalfTheta;  
72.     ratio2 = -0.333333333f + ratio2 * (sqU - 4.0f) * versHalfTheta;  
73.     ratio2 = 1.0f + ratio2 * (sqU - 1.0f) * versHalfTheta;  
74.   
75.     // Perform the bisection and resolve the folding done earlier.  
76.     f1 *= ratio1 * halfSecHalfTheta;  
77.     f2a *= ratio2;  
78.     f2b *= ratio2;  
79.     alpha *= f1 + f2a;  
80.     beta = f1 + f2b;  
81.   
82.     // Apply final coefficients to a and b as usual.  
83.     float w = alpha * q1w + beta * q2w;  
84.     float x = alpha * q1x + beta * q2x;  
85.     float y = alpha * q1y + beta * q2y;  
86.     float z = alpha * q1z + beta * q2z;  
87.   
88.     // This final adjustment to the quaternion’s length corrects for  
89.     // any small constraint error in the inputs q1 and q2 But as you  
90.     // can see, it comes at the cost of 9 additional multiplication  
91.     // operations. If this error-correcting feature is not required,  
92.     // the following code may be removed.  
93.     f1 = 1.5f - 0.5f * (w * w + x * x + y * y + z * z);  
94.     *dstw = w * f1;  
95.     *dstx = x * f1;  
96.     *dsty = y * f1;  
97.     *dstz = z * f1;  
98. }