1、旋轉的四元數表示
空間中的旋轉可用一個四元數來表述,如點 $P(\,x,y,z\,)$ 繞向量 $V(\,u,v,w\,)$ 旋轉 ,此旋轉過程可表示為四元數 Q $[\,\,cos( \frac{\theta}{2} ),\,sin( \frac{\theta}{2} )\cdot (u,v,w)\,]$ 。
2、基本旋轉的四元數表示
空間中的三維旋轉可視為繞三個基本軸的旋轉組合疊加,繞 $x,y,z$ 三個基本軸旋轉角度分別為 $ \phi,\theta,\psi$ ,則三個基本旋轉的四元素可表征為:
3、復合旋轉的四元數表示
繞三個基本軸的旋轉次序不同,其表征的空間旋轉也不同,下面給出順序為 $ZYX$ 、$XYZ$ 時的結果及相應推導過程。
3.1 對於 $ZYX$ ,其表征的旋轉過程定義為:
且已知四元數,可以反求歐拉角
3.2 對於 $XYZ$ ,其表征的旋轉過程定義為:
4、空間點旋轉表示
5、空間點旋轉的逆過程
對於坐標系中點 $P(\,x,y,z\,)$ ,其旋轉歐拉角為 $ (\,\phi,\theta,\psi\,)$ ,按照 $ZYX$ 順序旋轉后坐標點為 $P'(\,x',y',z'\,)$ ,完成 P 到 P’ 轉換過程,其逆過程,即 P’ 到 P 的變換則按照歐拉角 $ (\,-\phi,-\theta,-\psi\,)$,的順序進行。
6、程序實現
按照如上表述,編寫測試算例,展示一個橢圓顆粒的的空間變換旋轉過程。
1 clc;clear; 2 close all; 3 4 % 采用四元數方法定義旋轉 5 % 采用右手系,旋轉次序按照 Z-Y-X 進行 6 7 %% 原始橢圓形狀定義 8 psi = 0:pi/20:pi; 9 sita = 0:pi/20:2*pi; 10 a = 8; 11 b = 3.5; 12 c = 2; 13 x0 = 0; 14 y0 = 0; 15 z0 = 0; 16 17 x = zeros( length(psi), length(sita) ); 18 y = zeros( length(psi), length(sita) ); 19 z = zeros( length(psi), length(sita) ); 20 for j = 1 : length( psi ) 21 for i = 1 : length( sita ) 22 x(j,i) = x0 + a * sin( psi(j) ) * cos( sita(i) ); 23 y(j,i) = y0 + b * sin( psi(j) ) * sin( sita(i) ); 24 z(j,i) = z0 + c * cos( psi(j) ); 25 end 26 end 27 28 %% 橢圓旋轉 方法一:按照展開式,不借助函數, Z-Y-X 體系 29 sita_x = pi/3; 30 sita_y = pi/6; 31 sita_z = pi/5; 32 33 q = zeros(1,4); 34 q(1) = cos( 0.5 * sita_x ) * cos( 0.5 * sita_y ) * cos( 0.5 * sita_z ) + sin( 0.5 * sita_x ) * sin( 0.5 * sita_y ) * sin( 0.5 * sita_z ); 35 q(2) = sin( 0.5 * sita_x ) * cos( 0.5 * sita_y ) * cos( 0.5 * sita_z ) - cos( 0.5 * sita_x ) * sin( 0.5 * sita_y ) * sin( 0.5 * sita_z ); 36 q(3) = cos( 0.5 * sita_x ) * sin( 0.5 * sita_y ) * cos( 0.5 * sita_z ) + sin( 0.5 * sita_x ) * cos( 0.5 * sita_y ) * sin( 0.5 * sita_z ); 37 q(4) = cos( 0.5 * sita_x ) * cos( 0.5 * sita_y ) * sin( 0.5 * sita_z ) - sin( 0.5 * sita_x ) * sin( 0.5 * sita_y ) * cos( 0.5 * sita_z ); 38 39 T_rotate = zeros(3,3); 40 T_rotate(1,1) = q(1) * q(1) + q(2) * q(2) - q(3) * q(3) - q(4) * q(4); 41 T_rotate(1,2) = 2 * q(2) * q(3) - 2 * q(1) * q(4); 42 T_rotate(1,3) = 2 * q(1) * q(3) + 2 * q(2) * q(4); 43 T_rotate(2,1) = 2 * q(1) * q(4) + 2 * q(2) * q(3); 44 T_rotate(2,2) = q(1) * q(1) - q(2) * q(2) + q(3) * q(3) - q(4) * q(4); 45 T_rotate(2,3) = 2 * q(3) * q(4) - 2 * q(1) * q(2); 46 T_rotate(3,1) = 2 * q(2) * q(4) - 2 * q(1) * q(3); 47 T_rotate(3,2) = 2 * q(1) * q(2) + 2 * q(3) * q(4); 48 T_rotate(3,3) = q(1) * q(1) - q(2) * q(2) - q(3) * q(3) + q(4) * q(4); 49 50 x1 = zeros(size(x)); 51 y1 = zeros(size(y)); 52 z1 = zeros(size(z)); 53 54 for j = 1:size(x,1) 55 for i = 1:size(x,2) 56 p = [ x(j,i) - x0, y(j,i) - y0, z(j,i) - z0 ]; 57 x1(j,i) = x0 + T_rotate(1,1) * p(1) + T_rotate(1,2) * p(2) + T_rotate(1,3) * p(3); 58 y1(j,i) = y0 + T_rotate(2,1) * p(1) + T_rotate(2,2) * p(2) + T_rotate(2,3) * p(3); 59 z1(j,i) = z0 + T_rotate(3,1) * p(1) + T_rotate(3,2) * p(2) + T_rotate(3,3) * p(3); 60 end 61 end 62 63 %% 橢圓逆向歸位,方法一:按照展開式,不依賴函數,X-Y-Z 體系 64 sita_x = -sita_x; 65 sita_y = -sita_y; 66 sita_z = -sita_z; 67 68 q = zeros(1,4); 69 q(1) = cos( 0.5 * sita_x ) * cos( 0.5 * sita_y ) * cos( 0.5 * sita_z ) - sin( 0.5 * sita_x ) * sin( 0.5 * sita_y ) * sin( 0.5 * sita_z ); 70 q(2) = sin( 0.5 * sita_x ) * cos( 0.5 * sita_y ) * cos( 0.5 * sita_z ) + cos( 0.5 * sita_x ) * sin( 0.5 * sita_y ) * sin( 0.5 * sita_z ); 71 q(3) = cos( 0.5 * sita_x ) * sin( 0.5 * sita_y ) * cos( 0.5 * sita_z ) - sin( 0.5 * sita_x ) * cos( 0.5 * sita_y ) * sin( 0.5 * sita_z ); 72 q(4) = cos( 0.5 * sita_x ) * cos( 0.5 * sita_y ) * sin( 0.5 * sita_z ) + sin( 0.5 * sita_x ) * sin( 0.5 * sita_y ) * cos( 0.5 * sita_z ); 73 74 T_rotate = zeros(3,3); 75 T_rotate(1,1) = q(1) * q(1) + q(2) * q(2) - q(3) * q(3) - q(4) * q(4); 76 T_rotate(1,2) = 2 * q(2) * q(3) - 2 * q(1) * q(4); 77 T_rotate(1,3) = 2 * q(1) * q(3) + 2 * q(2) * q(4); 78 T_rotate(2,1) = 2 * q(1) * q(4) + 2 * q(2) * q(3); 79 T_rotate(2,2) = q(1) * q(1) - q(2) * q(2) + q(3) * q(3) - q(4) * q(4); 80 T_rotate(2,3) = 2 * q(3) * q(4) - 2 * q(1) * q(2); 81 T_rotate(3,1) = 2 * q(2) * q(4) - 2 * q(1) * q(3); 82 T_rotate(3,2) = 2 * q(1) * q(2) + 2 * q(3) * q(4); 83 T_rotate(3,3) = q(1) * q(1) - q(2) * q(2) - q(3) * q(3) + q(4) * q(4); 84 85 x3 = zeros(size(x)); 86 y3 = zeros(size(y)); 87 z3 = zeros(size(z)); 88 89 for j = 1:size(x,1) 90 for i = 1:size(x,2) 91 p = [ x1(j,i) - x0, y1(j,i) - y0, z1(j,i) - z0 ]; 92 x3(j,i) = x0 + T_rotate(1,1) * p(1) + T_rotate(1,2) * p(2) + T_rotate(1,3) * p(3); 93 y3(j,i) = y0 + T_rotate(2,1) * p(1) + T_rotate(2,2) * p(2) + T_rotate(2,3) * p(3); 94 z3(j,i) = z0 + T_rotate(3,1) * p(1) + T_rotate(3,2) * p(2) + T_rotate(3,3) * p(3); 95 end 96 end 97 98 %% 圖像顯示 99 figure 100 hold on 101 axis equal 102 axis off 103 104 quiver3( 0, 0, 0, 10, 0, 0, 'k', 'filled', 'linewidth', 3 ) 105 text( 10, 0, 0, '\it X', 'fontsize', 13, 'fontweight', 'bold' ) 106 quiver3( 0, 0, 0, 0, 12, 0, 'k', 'filled', 'linewidth', 3 ) 107 text( 0, 12, 0, '\it Y', 'fontsize', 13, 'fontweight', 'bold' ) 108 quiver3( 0, 0, 0, 0, 0, 10, 'k', 'filled', 'linewidth', 3 ) 109 text( 0, 0, 10, '\it Z', 'fontsize', 13, 'fontweight', 'bold' ) 110 111 % 原始橢圓 112 surf( x, y, z, 'facecolor', [0.8 0.8 0.8] ); 113 % 方法一旋轉 114 surf( x1, y1, z1, 'facecolor', [0.5 0.5 0.5] ); 115 % 方法一逆旋轉 116 surf( x3, y3, z3, 'facecolor', [0.7 0.5 0.3] );