[Insight-users] Transformation components from affine matrix

[Luis Ibanez]

Hi Hari,

You may want to look at

    Insight/Examples/Registration/ImageRegistration9.cxx:

lines 421-433:

   vnl_matrix<double> p(2, 2);
   p[0][0] = (double) finalParameters[0];
   p[0][1] = (double) finalParameters[1];
   p[1][0] = (double) finalParameters[2];
   p[1][1] = (double) finalParameters[3];
   vnl_svd<double> svd(p);
   vnl_matrix<double> r(2, 2);
   r = svd.U() * vnl_transpose(svd.V());
   double angle = asin(r[1][0]);

   std::cout << " Scale 1  = " << svd.W(0) << std::endl;
   std::cout << " Scale 2  = " << svd.W(1) << std::endl;
   std::cout << " Angle= " << angle * 45.0/atan(1.0) << std::endl;
	

Once you have the scaling and rotation components
of the matrix, you can compute the shearing as the
residual.


---

Another way of solving this is to use Polar Decomposition:
http://en.wikipedia.org/wiki/Polar_decomposition

This can be implemented as:

   //
   // 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".
   //
   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 )
       {
       // Temporary message here for debugging purposes...
       std::cout << "Polar decomposition used "
                 << ni << " iterations " << std::endl;
       break;
       }
     else
       {
       PQ = NQ;
       }
     }

   typename AffineTransformType::MatrixType QMatrix;

   QMatrix = NQ;

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



Where the QMatrix is the orthogonal component of the original matrix.



     Regards,


        Luis



-------------
Hari wrote:
> Hi all,
> 
>   Is there a way to get the different components (rotation, scaling and 
> shear) from a 2D affine matrix.
> I have a method for getting the rotation and scaling but am not sure how 
> to get the shear. The method I am using
> to get the rotation and scaling is:
> 
> Let A be the affine matrix.
> 
> USV' = svd(A);
> R = UV';
> 
> Scale X  = S[0,0]
> Scale Y = S[1,1]
> Angle = asin(R[1,0])
> 
> Is this the right way to do it ? Also I would appreciate if someone can 
> point me to how to get the shear factors.
> 
> Thanks!
> 
> 
> 
> 
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users