회전행렬에서 회전축'회전각의 분리(angular displacement)

http://www.devpia.com/Maeul/Contents/Detail.aspx?BoardID=30&MAEULNO=12&no=90&page=3

회전행렬에서 회전축'회전각의 분리(angular displacement)  | 3D Programming

2001-03-27 오전 3:46:24

이재권	번호: 90  / 읽음:3,999

안녕하세요.. 이 글은 회전행렬에서 회전축(고유값이 1인 고유벡터)과 그 축을 중심으로 한 회전각을 분리하는것에 대한 글입니다. 전에 쿼터니언을 꺼내 쓰다가 넘 오차가 많이 날때가 있어서 직접 풀어낸 공식입니다. 밤일(??) 하다 지쳐서(???) 여길 쫌 둘러봤는데.. 게임개발게시판은 상대적으로 글이 적네요^^: 그냥 혹시나 해서 올립니다. 일단 아래의 코드는 회전축과 회전각(쿼터니언 아님! 보통 angular displacement라고 도..)으로 행렬을 만드는 공식입니다. 이 공식을 유도하는 방법은 두가지죠. 보통 초보때 배우는 회전축->단위축일치->회전.... 이거하구 다른 하나는 쿼터니언-행렬 변환공식이죠. 둘 다 정리하면 아래의 공식이 나옵니다.   D3DMATRIX m;   float costheta = cos(theta), sintheta = sin(theta), Length;   // must be normalized if not   Length = sqrt(x*x + y*y + z*z);   x /= Length;   y /= Length;   z /= Length;   m._11 = (x*x) * (1-costheta) + costheta;   m._12 = (x*y) * (1-costheta) - (z*sintheta);   m._13 = (x*z) * (1-costheta) + (y*sintheta);   m._21 = (y*x) * (1-costheta) + (z*sintheta);   m._22 = (y*y) * (1-costheta) + costheta ;   m._23 = (y*z) * (1-costheta) - (x*sintheta);   m._31 = (z*x) * (1-costheta) - (y*sintheta);   m._32 = (z*y) * (1-costheta) + (x*sintheta);   m._33 = (z*z) * (1-costheta) + costheta; 보시면 알겠지만 행렬에 들어있는 값은 cos(theta/2), sin(theta/2)가 아닌 cos(theta), sin(theta)값이라서 여기서 곧바로 쿼터니언을 꺼낼 경우 대각원소의 절대값으로 기준행을 구한다 하더라도 오차가 쫌 나게 마련이드라구요.. 아래의 코드는 위의 공식을 역으로 풀어 정리한겁니다. 예전에 함 열라 연습장 더럽히면서 풀었던거죠^^: 아무리 찾아봐도 이에 관한 자료는 쿼터니언판밖에 없드라구여TT; void ADfromMAT(D3DMATRIX &m1, D3DVALUE &x, D3DVALUE &y, D3DVALUE &z, D3DVALUE  &theta) {   D3DVALUE _Trace, _sint, _cost, _cost1;   D3DVALUE   _11, _12, _13, _21, _22, _23, _31, _32, _33;     _11 = m1.m._11; _12 = m1.m._12; _13 = m1.m._13;     _21 = m1.m._21; _22 = m1.m._22; _23 = m1.m._23;     _31 = m1.m._31; _32 = m1.m._32; _33 = m1.m._33;     _Trace = _11 + _22 + _33;     if(_Trace == 3.0f) { // Trace = 1 + 2*cosθ 이므로 이 경우 θ가 0이 나온다.         x = 1.0f;         y = 0.0f;         z = 0.0f;         theta = 0.0f;         return;     }     _cost  = (_Trace - 1.0f)*0.5f;     if(_cost <= -1.0f) {         _cost = -1.0f;         theta = 3.141592654f;     } else if(_cost => 1.0f) {         _cost = 1.0f;         theta = 0.0f;     } else {         theta = acos(_cost);     } #ifdef AD_Fast_Extract        // ################################# method 1     if(_Trace+1.0f > 0.0f) { // |theta| < PIE         _Trace -= 1.0f;         _sint = sin(theta);  // theta값이 절대값이므로 결과값은 실제론 |sin (theta)|값이다.         x = (_32 - _23) / (_sint*2.0f);         y = (_13 - _31) / (_sint*2.0f);         z = (_21 - _12) / (_sint*2.0f);     } else { // |theta| = PIE         if( (_11 > _22) && (_11 > _33)) {  // use |x| for calc~             x = sqrt((_11 + 1.0f)*0.5f);             _Trace = 1.0f/(2.0f*x);             y = _12 * _Trace;             z = _13 * _Trace;         } else if(_22 > _33) {         // use |y| for calc~             y = sqrt((_22 + 1.0f)*0.5f);             _Trace = 1.0f/(2.0f*y);             x = _12 * _Trace;             z = _23 * _Trace;         } else {                   // use |z| for calc~             z = sqrt((_33 + 1.0f)*0.5f);             _Trace = 1.0f/(2.0f*z);             x = _13 * _Trace;             y = _23 * _Trace;         }         theta = ud_PI;     } #else   // ############################################### method 2         _cost1 = 1.0f / (1.0f - _cost);         if( (_11 > _22) && (_11 > _33)) {  // use |x| for calc~             x = sqrt((_11 - _cost)*_cost1);             _Trace = 1.0f/(2.0f*(1.0f - _cost)*x);             y = (_12 + _21) * _Trace;             z = (_13 + _31) * _Trace;             _Trace = 1.0f / (2.0f*x);             _sint = (_32 - _23)*_Trace;         } else if(_22 > _33) {         // use |y| for calc~             y = sqrt((_22 - _cost)*_cost1);             _Trace = 1.0f/(2.0f*(1.0f - _cost)*y);             x = (_12 + _21) * _Trace;             z = (_23 + _32) * _Trace;             _Trace = 1.0f / (2.0f*y);             _sint = (_13 - _31)*_Trace;         } else {                   // use |z| for calc~             z = sqrt((_33 - _cost)*_cost1);             _Trace = 1.0f/(2.0f*(1.0f - _cost)*z);             x = (_13 + _31) * _Trace;             y = (_23 + _32) * _Trace;             _Trace = 1.0f / (2.0f*z);             _sint = (_21 - _12)*_Trace;         }         if(_sint < 0.0f) theta = -theta; #endif } 우선 Trace를 이용하면 cosθ값이 간단하게 나온다. 그 다음 단위행렬이 아닐때 회전 각과 회전축을 구하는 방법은 두가지가 있다. ** method 1     theta값을 |theta| 로 고정시키고 x,y,z에 대한 연립방정식으로 푸는 방식이다.     계산속도는 빠르지만 회전각의 절대값이 π에 수렴할수록 회전축의 좌표값에서 오 차가 많이     난다. 회전각은 0 ~ 180도 범위로 나오며 회전축에 회전방향이 반영되어있다. ** method 2     x,y,z,sin(theta)에 대한 연립방정식으로 풀고 sin(theta)의 부호로 theta의 부호 를 결정한다.     이 경우엔 method 1보다 속도는 약간 느리지만 거의 오차가 발생하지 않는다.     회전각은 -180 ~ 180도 범위로 나온다. 공식유도가 잘 안될땐 한가지만 떠올려보세요.. 같은 회전을 의미하는 angular  displacement는 항상 2개라는걸..(단위행렬이 아니라면..) 사실 두가지 방법에서 속도차이는 거의 없습니다. 하지만 오차를 보면 method2가 훨씬 효율적이죠. 잘 알려진 보편적인 행렬->쿼터니언 변환공식은 method1과 같은방식입니다. 오차확대조건도 거의 비슷하죠.. 글쎄.. 이게 가장 좋은 방법이라곤 생각되지 않네요^^;  더 좋은 방법을 누군가 알고 계시다면 쫌 알려주시면 고맙겠습니다. 문의는 jacklee666@hotmail.com 으로 해 주시구여^^; 이 게시판에 쓸만한 정보가 많았으면 하네요.. 음.. 실은 요즘 2000용 서버작업중인데^^: 켕.. 여기서 CPIO에 관한 자료를 봐서리..  넘 고맙구 그냥 가기도 미안해서 이 글을 적습니다. 누구에겐가 조금이라도 도움이 되길 바랍니다. 수고하세요..