복소수의 $a + bi$형태처럼 사원수도 실수부와 허수부를 덧셈으로 묶고 실수를 사용해 허수부의 계수를 지정한다.
$$ q = a + bi + cj + dk $$
세 개의 허수부 $i$, $j$, $k$
복소수의 허수 단위 $i$와 같은 성질을 가진다.
$$ i^2 = j^2 = k^2 = -1 $$
두 허수의 곱은 나머지 하나의 허수에 대응 되는데, Chapter 10에서 본 회전 순환 순서와 유사하게 $i$ → $j$ → $k$ → $i$ 순서로 대응된다.
$$ ij = k $$
$$ jk = i $$
$$ ki = j $$
$$ ijk = jki = kij = -1 $$
Chapter 15에서 보았듯, 크기가 $1$이고 제곱한 수가 $-1$이 되는 허수 $i$와 같은 수가 있다면, 그것을 사용해 오일러 공식을 만족하게 할 수 있었다.
$$ e^{i\theta} = \cos\theta + i \sin\theta = (\cos\theta, \sin\theta)$$
사원수에서도 이런 수가 있다면 좋겠다.
크기가 $1$인 순허수 사원수를 고려해보자.
$$ q = (0, \vec{n}) $$
크기가 $1$이므로 다음과 같아서,
$$ \sqrt{0 + |\vec{n}|^2} = 1$$
하지만 막상 회전 사원수 $q$를 순허수 사원수 $\vec{v}$에 곱해보면, 그 결과가 순허수 사원수가 아니다.
$$ \vec{v’} = q \cdot \vec{v} $$
$$ = (\cos\theta, \sin\theta \cdot \vec{n}) \cdot (0, \vec{v}) $$
$$ = (-\sin\theta(\vec{n} \cdot \vec{v}), \cos\theta\vec{v} + \sin\theta(\vec{n} \times \vec{v})) $$
참고로, 회전 사원수 $q$를 왼쪽에 배치해야 한다. (순서가 중요하다.)
실수부가 $-\sin\theta(\vec{n} \cdot \vec{v})$ 라서 회전의 결과는 4차원에 존재하겠다.
사원수 곱의 결과가 항상 순허수 사원수가 되는 특별한 수식을 발견해보자.
회전 시킬 벡터 $\vec{v}$를 수직과 수평 성분으로 나누자.
$$\vec{v} = \vec{v_{\parallel}} + \vec{v_{\perp}}$$
이 벡터를 $\vec{n}$을 기준으로 각 $\theta$만큼 회전시킨 벡터는 어떻게 표현할 수 있을까?
수직인 성분만 회전되고, 수평인 성분은 그대로 두면 된다.
$$\vec{v’} = \vec{v_{\parallel}} + q \cdot \vec{v_{\perp}}$$
여기에 회전 사원수 $q$를 자연지수함수로 표현하면 다음과 같다.
$$\vec{v’} = \vec{v_{\parallel}} + e^{\vec{n}\theta} \cdot \vec{v_{\perp}}$$
(우변의 두번째 항) 회전 사원수와 직교하는 벡터의 곱의 경우
$$e^{\vec{n}\theta} \cdot \vec{v_{\perp}} = (0, \cos\theta\vec{v_{\perp}} + \sin\theta(\vec{n} \times \vec{v_{\parallel}}))$$
실수부가 $0$인 순허수 사원수가 되므로 3차원 공간의 벡터에 대응된다.
여기에 연산순서를 바꾼 후, 회전각에 $-\theta$를 대입하면 똑같은 결과가 나온다.
$$ \vec{v_{\perp}} \cdot e^{\vec{n}\theta} = (0, \cos\theta\vec{v_{\perp}} - \sin\theta(\vec{n} \times \vec{v_{\parallel}}))$$
$$ \vec{v_{\perp}} \cdot e^{\vec{n}(-\theta)} = (0, \cos\theta\vec{v_{\perp}} + \sin\theta(\vec{n} \times \vec{v_{\parallel}}))$$
따라서 직교하는 벡터의 다음과 같은 성질을 도출할 수 있다.
$$ e^{\vec{n}\theta} \cdot \vec{v_{\perp}} = \vec{v_{\perp}} \cdot e^{\vec{n}(-\theta)} $$
(수식에는 없었지만) 회전 사원수와 평행한 벡터의 곱의 경우
$$ e^{\vec{n}\theta} \cdot \vec{v_{\parallel}} = (-\sin\theta(\vec{n}\cdot \vec{v_{\parallel}}, \cos\theta\vec{v_{\parallel}})$$
사원수의 곱에서는 외적 계산 부분이 있었기 때문에 교환법칙이 성립하지 않았다. 그런데 이 경우 외적이 사라져서 교환법칙이 성립하게 된다.
$$ e^{\vec{n}\theta} \cdot \vec{v_{\parallel}} = \vec{v_{\parallel}} \cdot e^{\vec{n}\theta} $$
이러한 평행, 수식 벡터의 성질을 이용해서 이전에 보았던 회전을 전개하면 다음과 같이 풀어진다. (풀이 생략)
$$\vec{v’} = \vec{v_{\parallel}} + e^{\vec{n}\theta} \cdot \vec{v_{\perp}}$$
$$ = e^{\vec{n}(\frac{\theta}{2})} \cdot \vec{v} \cdot e^{\vec{n}(- \frac{\theta}{2})}$$
자연지수함수를 다시 사원수의 형태로 바꿔보자.
이번에는 회전 사원수 $q$를 $\theta$가 아닌 $\frac{\theta}{2}$만큼 회전하는 것으로 정의하자.
$$ q = (\cos\frac{\theta}{2}, \sin\frac{\theta}{2}\vec{n}) = e^{\vec{n}(\frac{\theta}{2})}$$
반대방향으로 회전하는 사원수는 켤레 사원수 이므로, 다음과 같이 최종 정리된다.
$$ e^{\vec{n}(\frac{\theta}{2})} \cdot \vec{v} \cdot e^{\vec{n}(- \frac{\theta}{2})} = q \cdot \vec{v} \cdot q^*$$
$q$를 간단히 다음과 같이 표현하자.
$$ q = (\cos\frac{\theta}{2}, \sin\frac{\theta}{2}\vec{n}) = (w, \vec{r})$$
여기서 $2(\vec{r} \times \vec{v})$를 $\vec{t}$로 치환하면 다음과 같이 간단하게 정리된다.
$$ \vec{v’} = \vec{v} + w\vec{t} + \vec{r} \times \vec{t} $$
만약 연속해서 다른 각도로 두 번 회전 하면 어떻게 될까?
첫 번째 각의 크기를 $2\alpha$, 두 번째 각의 크기를 $2\beta$로 하면, 연속으로 회전 변환하는 식은 다음과 같을 것이다.
$$ \vec{v’} = q_{\alpha} \cdot \vec{v} \cdot q_{\alpha}^* $$
$$ \vec{v’’} = q_{\beta} \cdot \vec{v’} \cdot q_{\beta}^* $$
두 번째 식에 첫 번째 식을 대입하면 다음과 같이 전개된다. (풀이 생략)
$$ \vec{v’’} = q_{\beta} \cdot (q_{\alpha} \cdot \vec{v} \cdot q_{\alpha}^* ) \cdot q_{\beta}^* $$
$$ = q_{(\alpha + \beta)} \cdot \vec{v} \cdot q_{(\alpha + \beta)}^* $$
따라서 두 차례 회전 사원수를 곱한 값은 두 각을 합한 회전 사원수의 값과 동일함을 알 수 있다.
예를 들어, 3차원 공간에서 $x$축을 중심으로 각 $\theta$만큼 회전시키는 데 사용하는 회전 사원수 $q$는 다음과 같을 것이다.
$$ q = (\cos\frac{\theta}{2}, \sin\frac{\theta}{2}\vec{x}) = \cos\frac{\theta}{2} + \sin\frac{\theta}{2}i $$
그러므로, 오일러 각을 구성하는 롤, 피치, 요 회전에 대응하는 각 사원수는 다음과 같이 표기할 수 있겠다.
$$ q_{roll} = \cos\frac{\theta_{roll}}{2} + \sin\frac{\theta_{roll}}{2}k $$
$$ q_{pitch} = \cos\frac{\theta_{pitch}}{2} + \sin\frac{\theta_{pitch}}{2}i $$
$$ q_{yaw} = \cos\frac{\theta_{yaw}}{2} + \sin\frac{\theta_{yaw}}{2}j $$
FORCEINLINEconstexprvoidQuaternion::FromRotator(constRotator&InRotator){floatsr=0.f,sp=0.f,sy=0.f;floatcr=0.f,cp=0.f,cy=0.f;// 오일러 각을 절반으로 줄이고 롤, 피치, 요의 sin, cos 값을 구한다.
Math::GetSinCos(sr,cr,InRotator.Roll*0.5f);Math::GetSinCos(sp,cp,InRotator.Pitch*0.5f);Math::GetSinCos(sy,cy,InRotator.Yaw*0.5f);// 위에서 살펴본 식에 따라서 W, X, Y, Z를 구한다.
W=cr*cp*cy+sr*sp*sy;X=cr*sp*cy+sr*cp*sy;Y=cr*cp*sy-sr*sp*cy;Z=sr*cp*cy-cr*sp*sy;}
위의 식에서 실수부를 $w$, 허수부의 값을 $x$, $y$, $z$로 치환한다. 요, 롤, 피치의 절반각은 $y$, $r$, $p$로 표기하고, 이들의 $\sin$값을 $sy$, $sr$, $sp$로 지정하고, $\cos$값을 $cy$, $cr$, $cp$로 지정하자.
$$ w = cr \cdot cp \cdot cy + sr \cdot sp \cdot sy $$
$$ x = cr \cdot sp \cdot cy + sr \cdot cp \cdot sy $$
$$ y = cr \cdot cp \cdot sy - sr \cdot sp \cdot cy $$
$$ z = sr \cdot cp \cdot cy - cr \cdot sp \cdot sy $$
롤 회전 값을 구해보자.
$2(wz + xy)$라는 값을 구해보면 다음과 같은 결과가 나온다. (풀이 생략)
$$ 2(wz + xy) = \sin\theta_{roll} \cdot \cos\theta_{pitch} $$
이번에는 $1 - 2(z^2 + x^2)$를 구해보면 다음과 같은 결과가 나온다.
$$ 1 - 2(z^2 + x^2) = \cos\theta_{roll} \cdot \cos\theta_{pitch} $$
여기서 $\sin\theta_{roll} \cdot \cos\theta_{pitch}$을 $\cos\theta_{roll} \cdot \cos\theta_{pitch}$로 나누면 $\tan\theta_{roll}$을 구할 수 있다.
$$ \frac{\sin\theta_{roll} \cdot \cos\theta_{pitch}}{\cos\theta_{roll} \cdot \cos\theta_{pitch}} = \tan\theta_{roll} = \frac{2(wz + xy)}{(1 - 2(z^2 + x^2)}$$
이제 $\arctan$ 함수로 롤 회전 값을 구할 수 있겠다.
$$ \theta_{roll} = atan2(2(wz + xy), 1 - 2(z^2 + x^2)) $$
피치 회전 값을 구해보자.
$wx - yx$라는 값을 구해보면 다음과 같은 결과가 나온다.
$$ wx - yx = \sin\frac{\theta_{pitch}}{2} \cdot \cos\frac{\theta_{pitch}}{2} $$
$$ = \frac{1}{2}\sin\theta_{pitch} $$
이제 $\arcsin$ 함수로 피치 회전 값을 구할 수 있겠다.
$$ \theta_{pitch} = asin(2(wx - yz)) $$
주의할 점은, Chapter 4에서 보았듯이 $\arcsin$ 함수는 정의역이 $[-1, 1]$ 범위로 제한되어 있으므로 $wx - yz$의 값이 $[-0.5, 0.5]$ 범위를 벗어나지 않는지 확인하는 것이 필요하다.
만약 벗어난다면, 범위 내 가장 가까운 값으로 설정해야 한다.
요 회전 값을 구해보자.
롤 회전처럼, 요 회전의 $\tan$ 함수 값을 구한 후 $\arctan$ 함수를 사용해서 구한다.
$$ \theta_{yaw} = atan2(2(wy + xz), 1 - 2(x^2 + y^2)) $$
구현 코드
FORCEINLINERotatorQuaternion::ToRotator()const{Rotatorresult;// (1) 2(wz + xy)
floatsinrCosp=2*(W*Z+X*Y);// (2) 1 - 2(z^2 + x^2)
floatcosrCosp=1-2*(Z*Z+X*X);// roll = atan2((1), (2))
result.Roll=Math::Rad2Deg(atan2f(sinrCosp,cosrCosp));// (3) wx - yz
floatpitchTest=W*X-Y*Z;// wx - yz은 [-0.5, 0.5]를 벗어나면 안 된다.
floatasinThreshold=0.4999995f;// 벗어나는 경우 범위 내 가장 가까운 값으로 설정한다.
if(pitchTest<-asinThreshold){result.Pitch=-90.f;}elseif(pitchTest>asinThreshold){result.Pitch=90.f;}// 벗어나지 않는 경우
else{// pitch = asin(2 * (3))
result.Pitch=Math::Rad2Deg(asinf(2*pitchTest));}// (4) 2(wy + xz)
floatsinyCosp=2*(W*Y+X*Z);// (5) 1 - 2(x^2 + y^2)
floatcosyCosp=1.f-2*(X*X+Y*Y);// yaw = atan2((4), (5))
result.Yaw=Math::Rad2Deg(atan2f(sinyCosp,cosyCosp));returnresult;}
위 행렬에서 트레이스 $t$를 구해서 $w$값을 알아내보자.
$$ t = 3 - 4(x^2 + y^2 + z^2) $$
여기에 1을 더하면 제곱근을 씌울 수가 있겠다.
$$ t + 1 = 4 - 4(x^2 + y^2 + z^2)$$
회전 사원수는 크기가 $1$이므로, $x^2 + y^2 + z^2$의 값은 $1 - w^2$가 된다.
$$ t + 1 = 4 - 4(1 - w^2) = 4w^2$$
이제 제곱근을 씌우자.
$$ \sqrt{t + 1} = 2w $$
따라서 사원수의 $w$값은 다음과 같이 구할 수가 있다.
$$ w = \frac{\sqrt{t + 1}}{2} $$
이제 $w$값을 이용해서 나머지 $x$, $y$, $z$값을 구해보자.
회전 변환행렬을 살펴보면, 대각 방향으로 대칭된 요소를 서로 빼면 두 사원수 요소의 곱으로 정리된다. 이것을 이용하자.
$$ 2yz + 2xw - 2yz + 2xw = 4xw $$
$$ 2xy + 2zw - 2xy + 2zw = 4zw $$
$$ 2xy + 2zw - 2xy + 2zw = 4zw $$
예외 사항
트레이스 값이 $t \leq -1$이면 제곱근을 구하는 값이 $0$보다 작거나 같아져서 $\sqrt{t + 1}$의 해가 존재하지 않게 된다.
이런 경우에는 켄 슈메이크(Ken Shoemake)의 알고리즘을 이용한다.
먼저 $w$를 제외하고, $x$, $y$, $z$ 중에서 가장 큰 요소를 파악하고 그 값을 구한다.
첫 번째 대각 요소, 두 번째 대각 요소, 세 번째 대각 요소를 서로 비교하면, 어느 것이 큰지 알 수 있다.
각 요소를 구하는 방법은 다음과 같다.
그다음 $xy$, $yz$, $xz$를 통해서 $x$, $y$, $z$의 값을 모두 구한다.
마지막으로 $w$값은 위에서 살펴봤던 $xw$, $yw$, $zw$를 통해서 구한다.
다음은 켄 슈메이크 알고리즘을 구현한 로직이다.
FORCEINLINEvoidQuaternion::FromMatrix(constMatrix3x3&InMatrix){floatroot=0.f;floattrace=InMatrix[0][0]+InMatrix[1][1]+InMatrix[2][2];if(trace>0.f){// W를 구한다.
root=sqrtf(trace+1.f);W=0.5f*root;// 나머지 X,Y,Z를 계산할 때 4W를 나눠주어야 한다.
root=0.5f/root;// X, Y, Z를 구한다.
X=(InMatrix[1][2]-InMatrix[2][1])*root;Y=(InMatrix[2][0]-InMatrix[0][2])*root;Z=(InMatrix[0][1]-InMatrix[1][0])*root;}else{BYTEi=0;// X,Y,Z 중에서 가장 큰 요소를 파악해서 행렬상의 위치를 i에 담는다.
if(InMatrix[1][1]>InMatrix[0][0]){i=1;}if(InMatrix[2][2]>InMatrix[i][i]){i=2;}// i, j, k 의 순서를 지정한다.
// x, y, z를 구할 때 필요한 행렬상의 위치가 된다.
staticconstBYTEnext[3]={1,2,0};BYTEj=next[i];BYTEk=next[j];// 가장 큰 요소의 값을 i, j, k 를 이용해서 구한다.
root=sqrtf(InMatrix[i][i]-InMatrix[j][j]-InMatrix[k][k]+1.f);float*qt[3]={&X,&Y,&Z};// 가장 큰 요소를 저장한다.
*qt[i]=0.5f*root;root=0.5f/root;// 나머지 두 요소의 값을 구해서 저장한다.
*qt[j]=(InMatrix[i][j]+InMatrix[j][i])*root;*qt[k]=(InMatrix[i][k]+InMatrix[k][i])*root;// 마지막으로 W 값을 구한다.
W=(InMatrix[j][k]-InMatrix[k][j])*root;}}
선형 보간은 빠르고 간편하나, 원의 궤적을 따라 발생하는 회전을 정확히 반영하지는 못한다.
두 사원수가 이루는 각을 $\theta$라고 할 때, 주어진 비율 $t$에 따라 각각 $t\theta$, $(1-t)\theta$로 나눠서 구할 수 있겠다.
이것을 만족시키는 두 계수 $\alpha$와 $\beta$를 찾아야 한다.
$$ q’ = \alpha \cdot q_1 + \beta \cdot q_2 $$
우리가 구하려는 벡터를 $\vec{v}$는 서로 직교하는 두 기저벡터 $\vec{x}$와 $\vec{y}$의 조합으로 다음과 같이 계산할 수 있다.
$$ \vec{v} = \cos(t\theta)\vec{x} + \sin(t\theta)\vec{y} $$
여기서 우리가 모르는 $\vec{y}$를 구해야 하겠다.
위 그림을 보면, $\vec{u} - \cos\theta\vec{x}$와 $\sin\theta\vec{y}$의 값은 같다. 따라서 다음과 같이 $\vec{y}$를 구할 수 있겠다.
$$ \vec{y} = \frac{\vec{u} - \cos\theta\vec{x}}{\sin\theta} $$
이것을 $\vec{v}$ 식에 대입하면 다음과 같이 전개 된다.
$$ \vec{v} = \cos(t\theta)\vec{x} + \sin(t\theta)(\frac{\vec{u} - \cos\theta\vec{x}}{\sin\theta}) $$
$$ = \frac{\sin((1 - t)\theta)}{\sin\theta}\vec{x} + \frac{\sin(t\theta)}{\sin\theta}\vec{y} $$
따라서 이렇게 구한 두 계수 $\alpha$와 $\beta$를 $ q’ = \alpha \cdot q_1 + \beta \cdot q_2 $ 식에 대입하면 최종 식이 만들어진다.
$$ q’ = \frac{\sin((1 - t)\theta)}{\sin\theta} \cdot q_1 + \frac{\sin(t\theta)}{\sin\theta} \cdot q_2 $$
고려할 점 (1)
두 벡터가 만들어내는 회전 경로는 언제나 짧은 경로와 긴 경로 두 가지가 존재한다.
만약에 언제나 짧은 경로를 사용하고 싶다면?
두 벡터의 사잇각이 $180\degree$보다 작은지 확인해야 할 것이다.
하지만 우리는 사원수를 사용하고 있고, 사원수는 각도의 절반을 저장하고 있다. 그렇기 때문에 예를 들어, 저장된 값이 $120\degree$ 라면 회전 값은 그 두 배인 $240\degree$가 된다.
따라서 짧은 경로인지 판단하는 값은 $90\degree$가 되어야 하겠다.
두 사원수의. 내적값이 음수라면 사잇각이 $90\degree$보다 크다는 의미이므로, 이 때 짧은 경로의 사원수로 바꿔주면 된다.
짧은 경로의 사원수는 사원수를 반전시키는 것으로 간단히 구할 수 있다.
예를 들어, 사잇각이 $120\degree$인 경우에 실제 회전은 $240\degree$가 된다. 이 회전과 똑같아지는 짧은 경로의 회전은 $-120\degree$ 회전이다. 그리고 $-120\degree$에 대응하는 사원수는 그 절반인 $-60\degree$을 가진 사원수가 된다. 이것은 해당 사원수를 반전 시켜서 얻을 수 있다.
고려할 점 (2)
두 사원수의 방향이 평행하면 분모의 $\sin\theta$의 값이 $0$이 되어서 계산이 불가능해진다.
따라서 이 경우에는 선형 보간을 하도록하자.
구현 코드
FORCEINLINEQuaternionQuaternion::Slerp(constQuaternion&InQuaternion1,constQuaternion&InQuaternion2,floatInRatio){// 입력된 두 개의 사원수 q1, q2
Quaternionq1=InQuaternion1,q2=InQuaternion2;// 두 사원수의 내적을 구한다.
floatdot=q1.X*q2.X+q1.Y*q2.Y+q1.Z*q2.Z+q1.W*q2.W;// 내적 값이 0보다 작으면 짧은 거리를 사용하도록 방향을 전환한다.
if(dot<0.0f){q1=-q1;dot=-dot;}floatalpha=1.f,beta=0.f;// 두 사원수의 사잇각이 작으면 선형 보간으로 수행한다.
if(dot>0.9995f){alpha=1.0f-InRatio;beta=InRatio;}else{// 내적값의 아크코사인으로 theta를 구한다.
constfloattheta=acosf(dot);// 분모를 계산한다.
constfloatinvSin=1.f/sinf(theta);// alpha, beta를 구한다.
alpha=sinf((1.f-InRatio)*theta)*invSin;beta=sinf(InRatio*theta)*invSin;}// 최종 사원수를 구한다.
Quaternionresult;result.X=alpha*q1.X+beta*q2.X;result.Y=alpha*q1.Y+beta*q2.Y;result.Z=alpha*q1.Z+beta*q2.Z;result.W=alpha*q1.W+beta*q2.W;returnresult;}
