[Insight-users] Problems with VersorRigid3Transform
Luis Ibanez
luis.ibanez@kitware.com
Wed May 12 15:25:18 EDT 2004
Hi Anders,
You are not taking advantage of the beauty of Versors :-)
Describing rotations as a combinations of rotations around
the three axis is a *very bad* approach. E.g. having a
rotation that is 30 degrees around X, then 30 degrees around
Y, and finally 30 degrees around Z... this is only good for
playing with the Rubirk's cube. :-/
Since 1866 Hamilton came up with Versors for representing
rotations, and there were very good reasons for that.
You should think of a rotation as something defined by a
rotation-axis in 3D space (not the axis of coordinates),
and an angle. It is a single operation.
--
About the problems in your code:
1) You seem to be expecting that the Transform will
accumulate Translations. This is not done in ITK.
When you set the Translation of the VersorRigid3DTransform
you don't compose with the previous translation, you
replace the previous translation
2) The same goes about rotations. When you call SetRotation(),
you are not composing with previous rotations, you simply
replace any pre-existing rotation.
3) The Versor Rigid 3D transform do all the CenterOfRotation
calculations for you. Please remove all the image center
computation and simply call SetCenter(). You may find
interesting to look at the source code in
Insight/Code/Common/
itkVersorRigid3DTransform.h
itkVersorRigid3DTransform.txx
4) When you get the parameters array from this transform,
the first three values are the x,y,z components of the
Versor, while the last three values are the x,y,z components
of the translation. The x,y,z components of the versor
are equal to a unit vector multiplied by sin(angle/2).
Please read the tutorials on Quaternions in:
http://www.itk.org/cgi-bin/cvsweb.cgi/InsightDocuments/Developer/General/QuaternionsI.pdf?cvsroot=Insight
http://www.itk.org/cgi-bin/cvsweb.cgi/InsightDocuments/Developer/General/QuaternionsII.pdf?cvsroot=Insight
Note that in general you don't want to use Quaternions
but Versors. A Versor is the equivalent of a Unit Quaternion.
A full Quaternion could represent rotations and scaling, but
you are usually interested only in the rotational part, that
is what the Versor represent.
Regards,
Luis
--------------------
Anders Almås wrote:
> Hi!
>
> I am performing registration of 3D images and a I am using
> VersorRigid3DTransform. My goal is to create a test set by using
> VersorRigid3DTransform and then try to registrate it back to the
> starting point.
>
> My problem is to do the set up for this test... It might be a dumb
> question but...
>
> My code is the following:
>
>
> typedef itk::ResampleImageFilter<ImageType,ImageType> FilterType;
> //The FilterType is for checking the transform/rotation..
>
> FilterType::Pointer filter = FilterType::New();
>
> const int offsetx = par[0]; //Translation parameters in
> const int offsety = par[1]; // x-, y- and z direction
> const int offsetz = par[2];
>
>
> typedef itk::VersorRigid3DTransform<double> TransformType;
>
> TransformType::Pointer rotation = TransformType::New();
> TransformType::VersorType VersorType;
>
> typedef VersorType::VectorType VectorType;
>
> TransformType::Pointer transform = TransformType::New();
> TransformType::OutputVectorType translation;
> translation[0] = offsetx;
> translation[1] = offsety;
> translation[2] = offsetz;
>
> transform->SetTranslation(translation);
>
> typedef itk::LinearInterpolateImageFunction<ImageType,double>
> InterpolateType;
>
> InterpolateType::Pointer interpolator = InterpolateType::New();
>
> filter->SetInterpolator(interpolator);
>
> filter->SetDefaultPixelValue(0);
>
> double spacing[ImageDimension]; // 3D...
>
> spacing[0] = 1.453; // Spacings between the pixels..
> spacing[1] = 1.453;
> spacing[2] = 2.0;
>
> filter->SetOutputSpacing(spacing);
>
> double origin[ImageDimension];
>
> origin[0]=0;
> origin[1]=0;
> origin[2]=0;
>
> filter->SetOutputOrigin(origin);
>
> ImageType::SizeType size;
>
> size[0] = input->GetBufferedRegion().GetSize()[0];
> size[1] = input->GetBufferedRegion().GetSize()[1];
> size[2] = input->GetBufferedRegion().GetSize()[2];
>
> filter->SetSize(size);
>
>
> filter->SetInput(input);
> TransformType::OutputVectorType translation1;
>
> //Calculate a imagecenter..
>
> const double imageCenterX = origin[0] + spacing[0] * size[0] / 2.0;
> const double imageCenterY = origin[1] + spacing[1] * size[1] / 2.0;
> const double imageCenterZ = origin[2] + spacing[2] * size[2] / 2.0;
>
> translation1[0] = -imageCenterX;
> translation1[1] = -imageCenterY;
> translation1[2] = -imageCenterZ;
>
> transform->SetTranslation( translation1 );
>
> double pi = 3.14159265358979323846;
>
> const double degreesToRadians = pi/180;
> const double angle1 =par[3] * degreesToRadians; //par[3] =
> rotation in x direction
>
> VersorType rotation_x;
> VersorType rotation_y;
> VersorType rotation_z;
>
>
> VectorType axis_x;
> axis_x[0] = 1; // Rotate around x-axis...
> axis_x[1] = 0;// Rotate around y-axis...
> axis_x[2] = 0;// Rotate around z-axis...
>
> rotation_x.Set(axis_x,angle1);
> transform->SetRotation(rotation_x);
>
> const double angle2 =par[4] * degreesToRadians; //Rotation in
> y-direction
>
> VectorType axis_y;
> axis_y[0] = 0; // Rotate around x-axis...
> axis_y[1] = 1;// Rotate around y-axis...
> axis_y[2] = 0;// Rotate around z-axis...
>
> rotation_y.Set(axis_y,angle2);
> transform->SetRotation(rotation_y);
>
> const double angle3 =par[5] * degreesToRadians; // Rotation in z
> direction
> VectorType axis_z;
> axis_z[0] = 0; // Rotate around x-axis...
> axis_z[1] = 0;// Rotate around y-axis...
> axis_z[2] = 1;// Rotate around z-axis...
>
> rotation_z.Set(axis_z,angle3);
> transform->SetRotation(rotation_z);
>
> TransformType::OutputVectorType translation2;
> translation2[0] = imageCenterX;
> translation2[1] = imageCenterY;
> translation2[2] = imageCenterZ;
> transform->SetTranslation( translation2); //Translate back
>
> filter->SetTransform( transform );
> filter->Update();
> // TransformType::MatrixType matrix = transform->GetRotationMatrix();
> TransformType::ParametersType param = transform->GetParameters();
> std::cout<<"Transformen: "<<param<<std::endl;
> // return filter->GetOutput();
>
>
> What is wrong.. when i execute the code, it seems it is only performing
> rotation in the z - direction. Because when I execute the result with
> parameters: The rotation in all directions is 30 degrees.. and the
> translation is 10 mm in each direction...
> the result is:
> [0,0,0,0.258819,185.984,185.984,176]
>
> if I have understand it correctly:
>
> the three first zero's are versors?? and 0.258819 is the total
> rotation?? or what?? and the final three are the image center. Is this
> correct?
> Do I now have an image which is rotated in 30 degrees in every direction
> and translated 10 mm in the same directions?? My second question is:
> how to apply these parameters on the resampling filter?
>
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users@itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
More information about the Insight-users
mailing list