[Insight-users] Volume not perfectly rotated

[Luis Ibanez]

Hi Nic,


It seems that you have not initialized the Center of Rotation of the
VersorRigid3DTransform.

By default that center is at (0,0,0), so it is likely that your image
if being flipped around the corner of the image, not around the center
of the image.

You may want to compute the center of the image in physical coordinates
and use that as the center of rotation for the VersorRigid3DTransform.

You will find examples on how to compute the center of the image in
the ITK Software Guide.

       http://www.itk.org/ItkSoftwareGuide.pdf


It comes down to

     for( i=0; i<3; i++)
       {
       Center[i] = Origin[i] + Size[i] * Spacing[i] / 2.0
       }

-----


Also, Please aware that *you are not using the Quaternion correctly*.


   //: Construct quaternion from Euler Angles,
   // That is a rotation about the X axis, followed by Y, followed by
   // the Z axis, using a fixed reference frame.
   vnl_quaternion(T theta_X, T theta_Y, T theta_Z);



That is not the order of the rotational sequence that you are assuming.



It is too bad that so many people have grown used to Euler angles.
They are one of the *worst things* that have made their way in
computer graphics. They are clumsy and provide a very poor
representation of the SO(3) rotational space

          http://en.wikipedia.org/wiki/SO%283%29

Euler angles have impoverished the understanding of generations of
engineers on the nature of rotational space and crippled them to
think in Cartesian schems instead of embracing the natural properties
of the rotational space. That results in uncountable pieces of software
that require "if" conditions when computing rotational components.



---------

The correct way of composing rotations in space when using
Quaternions is to define a quaternion for each one of the
intermediate rotations and then use the Quaternion composition
method:

     vnl_quaternion::operator*()   in VNL

or

     itkVersor::operator*()   in ITK


In your case:

 > What I do is simply initialize the versorRigid3dTransform with a 180°
 > rotation around y-axis and -94.55 rotation around z-axis.


You should do:

typedef itk::Versor<double>      VersorType;
typedef VersorType::VectorType   AxisType;
typedef VersorType::ValueType    ValueType;

AxisType axis1;   // Rotation Axis = Y
axis1[0] = 0.0;
axis1[0] = 1.0;
axis1[0] = 0.0;

ValueType angle1 = 3.141516...

VersorType versor1;
versor1.Set( axis1, angle1 ); // 180 degrees around Y


AxisType axis2;   // Rotation Axis = Z
axis2[0] = 0.0;
axis2[0] = 0.0;
axis2[0] = 1.0;

ValueType angle2 = 3.141516 * ( -94.55 ) / 180.0;

VersorType versor2;
versor2.Set( axis2, angle2 ); // -94.55 degrees around Z


      VersorType versor3 = versor2 * versor1;


The last variable "versor3" will contain the correct rotation
expressed in terms of a unit quaternion.


Note that most of the time, when people use Quaternions,
what they actually want to use are Versors.

Quaternions can represent rotation *and scaling* in space,
while Versors are limited to rotations. Versors are equivalent
to Unit-Quaternions.


You may want to read more about Quaternions and their properties:
http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation

and the ITK Quaternion Tutorials:

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


    Regards,


       Luis



----------------
Nic wrote:
> Hi all,
>     I have a strange behaviour I would like to correct in my initialisation.
> I used quaternion for rotation initialisation like this, as previously
> adviced :
> 
> *******************************************************
> double RotX = angleX * vnl_math::pi / 180.0;
> double RotY = angleY * vnl_math::pi / 180.0;
> double RotZ = angleZ * vnl_math::pi / 180.0;
> vnl_quaternion<double> RR(RotX,RotY,RotZ);
> rotation.Set(RR);
> transform->SetRotation( rotation );
> 
> [...]
> // Checking
> VersorType rotationInitiale = transform->GetVersor();
> vnl_vector_fixed<double,3> anglesEuleriensInitiaux = 
> (rotationInitiale.GetVnlQuaternion()).rotation_euler_angles();
> cout << "Angles d'Euler (initiaux): " << 
> (anglesEuleriensInitiaux*180.0/vnl_math::pi) << endl;
> 
> ********************************************************
> 
> What I doesn't understand actually is why by flipping around X or Y axis 
> by 180°, I get an "half image" at the beginning of the stack
> Is this behavious linked to the Quaternions ? Is there a special case 
> for a pi angle ? Is there a way to avoid this behaviour ?
>  
> Images:
> http://itk.fete.ch/perso/Rot_0_0_0.jpg
> http://itk.fete.ch/perso/Rot_0_180_0.jpg
> http://itk.fete.ch/perso/Rot_180_0_0.jpg
>  
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users