[Insight-users] Get angle values from a 3D affine transform

[Luis Ibanez]

Hi Yan,

A better method for extracting rotation information is
to compute the Polar decomposition of the Affine Matrix.

For background in this method, please look at:

Polar decomposition algorithm proposed by [Higham 86]
SIAM J. Sci. Stat. Comput. Vol. 7, Num. 4, October 1986.
"Computing the Polar Decomposition - with Applications"
by Nicholas Higham.

recommended by
Shoemake in the paper "Matrix Animation and Polar Decomposition".
http://citeseer.ist.psu.edu/shoemake92matrix.html



You can do this with the following code:



   typedef vnl_matrix<double> VnlMatrixType;

   VnlMatrixType M = initialAffine->GetMatrix().GetVnlMatrix();

   VnlMatrixType PQ = M;
   VnlMatrixType NQ = M;
   VnlMatrixType PQNQDiff;

   const unsigned int maximumIterations = 100;

   for(unsigned int ni = 0; ni < maximumIterations; ni++ )
     {
     // Average current Qi with its inverse transpose
     NQ = ( PQ + vnl_inverse_transpose( PQ ) ) / 2.0;
     PQNQDiff = NQ - PQ;
     if( PQNQDiff.frobenius_norm() < 1e-7 )
       {
       std::cout << "Polar decomposition used "
                 << ni << " iterations " << std::endl;
       break;
       }
     else
       {
       PQ = NQ;
       }
     }

   MatrixType QMatrix;

   QMatrix = NQ;

   std::cout << "Initial Matrix = " << std::endl << M << std::endl;
   std::cout << "Q Matrix = " << std::endl << QMatrix << std::endl;

   VersorType versor;
   versor.Set( QMatrix );

   std::cout << versor.GetAxis() << std::endl;
   std::cout << versor.GetAngle() << std::endl;



Note that Versors (=unit Quaternions) are a much better description
of rotations in 3D space than Euler angles. Euler angles are ambiguous
unles you specify the order in which the rotations are applied.

Note also that the components of a versor *are not* equivalent to
Euler angles.  For details on Quaternions you may want to look at
the ITK Tutorials:

       http://www.itk.org/HTML/Tutorials.htm

in particular to:

  http://www.itk.org/CourseWare/Training/QuaternionsI.pdf
  http://www.itk.org/CourseWare/Training/QuaternionsII.pdf



  Regards,


     Luis




---------------
Yan Yang wrote:
> Hi All,
>  
> Could you have a look at my method, I am not sure if it is correct. 
> Thanks in advance.
>  
> I am trying to extract angle information from a 3D affine transform, I 
> want to know if I can do as the following steps.
>  
> 1. Decompose the affine transform matrix (A) excluding translation using 
> SVD (vnl_svd), get A = UWV,
>     so the rotaiton matrix can be R = UV'.
>  
> 2. Create an euler transform (E) using the rotation matrix, 
> E->setMatrix( R ),
>  
> 3. Calculate the three angles around axes by
>     angleX = E->GetAngleX()
>     angleY = E->GetAngleY()
>     angleZ = E->GetAngleZ()
>  
> Am I right?
>  
> By the way, should the rotation matrix R be vnl_matrix type? Does 
> vnl_matrix support 3D?
>  
> Thanks,
>  
> Yan
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users