[Insight-users] Volume not perfectly rotated
Nic
itk at fete.ch
Thu Oct 11 12:52:04 EDT 2007
Hi Luis,
thanks for those clarifications. I forgot some precisions :)
Firstly, I didn't paste the setcenter part of my code, that was above the pasted code.
As a remember, I have two stacks resulting from a "top-down scan" and a "bottom-up" scan, both from the outside of the sample until a little bit further than the sample middle.
Each stack contains an overlapping volume, identified by :
- firstFixedOverlappingBlockImage to lastFixedImage in "top stack"
- firstMovingOverlappingBlockImage to lastMovingImage in "bottom stack"
Thus I can easily geometrically calculate the center of the overlapping volume in the fixed stack and the same for the moving stack.
Translation is then set to (centerMovingImage-centerFixedImage) and rotation center set to the fixedOverlappingVolumeCenter.
Moreover, for translation fine tuning reasons, I allow an offsetX,Y,Z to be added when calculating the center of the moving image overlapping block. Here is the code..
********************************************************************************************
// ----------------------------------------
// OVERLAPPING BLOCKS CENTER COMPUTATION
// ----------------------------------------
// FIXED IMAGE
const SourceImageType::SpacingType fixedSpacing = pile3Drecto->GetSpacing();
double fixedImageOverlappingBlockCoordinates[3];
fixedImageOverlappingBlockCoordinates[0] = fixedImageSize[0] / 2.0;
fixedImageOverlappingBlockCoordinates[1] = fixedImageSize[1] / 2.0;
fixedImageOverlappingBlockCoordinates[2] = (lastFixedImage+firstFixedOverlappingBlockImage)/ 2.0;
TransformType::InputPointType centerFixedImage;
centerFixedImage[0] = fixedOrigin[0]+ fixedSpacing[0] * fixedImageOverlappingBlockCoordinates[0];
centerFixedImage[1] = fixedOrigin[1]+ fixedSpacing[1] * fixedImageOverlappingBlockCoordinates[1];
centerFixedImage[2] = fixedOrigin[2]+ fixedSpacing[2] * fixedImageOverlappingBlockCoordinates[2];
// MOVING IMAGE
const SourceImageType::SpacingType movingSpacing = pile3Dverso->GetSpacing();
double movingImageOverlappingBlockCoordinates[3];
movingImageOverlappingBlockCoordinates[0] = (movingImageSize[0] / 2.0)+offsetX;
movingImageOverlappingBlockCoordinates[1] = (movingImageSize[1] / 2.0)+offsetY;
movingImageOverlappingBlockCoordinates[2] = ((lastMovingImage+firstMovingOverlappingBlockImage)/ 2.0)+offsetZ;
TransformType::InputPointType centerMovingImage;
centerMovingImage[0] = movingOrigin[0]+ movingSpacing[0] * movingImageOverlappingBlockCoordinates[0];
centerMovingImage[1] = movingOrigin[1]+ movingSpacing[1] * movingImageOverlappingBlockCoordinates[1];
centerMovingImage[2] = movingOrigin[2]+ movingSpacing[2] * movingImageOverlappingBlockCoordinates[2];
cout << "Centre de rotation MOVING (coordonnees): [" << movingImageOverlappingBlockCoordinates[0] << ";" <<movingImageOverlappingBlockCoordinates[1] << ";" << movingImageOverlappingBlockCoordinates[2] << "]" << endl;
cout << "Centre de rotation MOVING (um): [" << centerMovingImage[0] << ";" << centerMovingImage[1] << ";" << centerMovingImage[2] << "]" << endl;
cout << "OFFSET (coordonnees): [" << offsetX << ";" << offsetY << ";" << offsetZ << "]" << endl;
// -------------------------------------------------------------------------
// SET ROTATION CENTER SET AND TRANSLATION VECTOR
// -------------------------------------------------------------------------
// set center of fixed image as the rotation center
transform->SetCenter( centerFixedImage );
// translation set as the vector relating the center of the moving image to the center of the fixed image
transform->SetTranslation( centerMovingImage-centerFixedImage);
********************************************************************************************
For the quaternion computation, in fact my real code is a little bit different than the bad test code I sent.
As you said, I constructed the quaternion for each rotation, multiplied and set the rotation like this:
*************************************************************
typedef TransformType::VersorType VersorType;
VersorType rotation;
// rotation d'un angle angleZ selon Oz, dans les deux cas de flip
double RotZ = angleZ * vnl_math::pi / 180.0;
vnl_quaternion<double> Rz(0, 0, RotZ);
vnl_quaternion<double> RFlip, RR;
if (flipHorizVert) RFlip = vnl_quaternion<double>(0, vnl_math::pi, 0); // rotation de 180° selon Oy
else RFlip = vnl_quaternion<double>(vnl_math::pi, 0, 0); // rotation de 180° selon Ox
RR = RFlip * Rz;
rotation.Set(RR);
transform->SetRotation( rotation );
registration->SetInitialTransformParameters( transform->GetParameters() );
*************************************************************
All this code gave a quite nice initialization result:
- 3D view: http://itk.fete.ch/perso/3DViewInit.jpg
- Slice view: http://itk.fete.ch/perso/SliceInit001.jpg and http://itk.fete.ch/perso/SliceInit002.jpg
But as you can see on http://itk.fete.ch/perso/SliceInit003.jpg, the first (red) image of the moving stack via the initialized transform is half cut.
I'll change my code for your versor code and see if it happend again.
Cheers, nic
----- Original Message -----
From: "Luis Ibanez" <luis.ibanez at kitware.com>
To: "Nic" <itk at fete.ch>
Cc: "Insight Users" <insight-users at itk.org>
Sent: Thursday, October 11, 2007 5:44 PM
Subject: Re: [Insight-users] Volume not perfectly rotated
>
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20071011/f0a64a7f/attachment-0001.html
More information about the Insight-users
mailing list