[Insight-users] Affine transform based registration in 3D

[Luis Ibanez]

Hi Marco,


You got most of it right.

The rotation aspect, however, requires a bit more work.

The AffineTransform can give you a combination of

* Shearing
* Scaling
* Rotation

So, first you must separate the rotational part of the Matrix
from the rest (scaling and shearing). This can be done by
using Polar decomposition,  for Details, please See:
http://www.cmake.org/pipermail/insight-users/2006-August/019025.html

Once you extract the Orthogonal matrix (that represents
rotation) out of the Affine Matrix, you can feed it into a
Versor, by using the method:

                 versor.SetMatrix( orthogonalMatrix )

and finally, from the Versor, you can get:

        double angle = versor.GetAngle();
        VectorType rotationAxis = versor.GetAxis();

that will tell you how many *radians* is the current Affine
transform rotating, and around what axis in space.


      Regards,


            Luis


-----------------------------------
On Fri, Dec 11, 2009 at 9:41 AM, marco giordano <marco.giord at gmail.com> wrote:
> Hi,
>
> I am running an Affine transform registration in 3D (basically
> ImageRegistration9.cxx).
>
> In order to support 3D mode, appart from some obvious changes, I did the
> following:
>
> Extend the optimizer scale:
>
> ....
>  optimizerScales[0] =  1.0;
>   optimizerScales[1] =  1.0;
>   optimizerScales[2] =  1.0;
>   optimizerScales[3] =  1.0;
>   optimizerScales[4] =  1.0;
>   optimizerScales[5] =  1.0;
>   optimizerScales[6] =  1.0;
>   optimizerScales[7] =  1.0;
>   optimizerScales[8] =  1.0;
>   optimizerScales[9] =  translationScale;
>   optimizerScales[10] =  translationScale;
>   optimizerScales[11] =  translationScale;
> ...
>
>
> Get the full set of final parameters to print out
>
> ....
> const double finalRotationCenterX = transform->GetCenter()[0];
>   const double finalRotationCenterY = transform->GetCenter()[1];
>   const double finalRotationCenterZ = transform->GetCenter()[2];
>   const double finalTranslationX    = finalParameters[9];
>   const double finalTranslationY    = finalParameters[10];
>   const double finalTranslationZ    = finalParameters[11];
> ....
>
> Then I guess I must retrieve the rotation angle by SVD
> which I guess it is run on the 3x3 matrix but I do not know in which order I
> must pick up the value. I am not sure that I got
> this right but my guess would be:
>
> (Guess)
>  vnl_matrix<double> p(3, 3);
>   p[0][0] = (double) finalParameters[0];
>   p[0][1] = (double) finalParameters[1];
>   p[0][2] = (double) finalParameters[2];
>
>   p[1][0] = (double) finalParameters[3];
>   p[1][1] = (double) finalParameters[4];
>   p[1][2] = (double) finalParameters[5];
>
>   p[2][0] = (double) finalParameters[6];
>   p[2][1] = (double) finalParameters[7];
>   p[2][2] = (double) finalParameters[8];
>
>
>   vnl_svd<double> svd(p);
>   vnl_matrix<double> r(3, 3);
>   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 << " Scale 3         = " << svd.W(3)                 <<
> std::endl;
>   std::cout << " Angle (degrees) = " << angle * 45.0 / atan(1.0) <<
> std::endl;
>
>
> Does that make sense ?
>
> I would appreciate some help on that.
>
> Thank you
>
> --
> Marco G.
>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
>