<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=Content-Type content="text/html; charset=iso-8859-1">
<META content="MSHTML 6.00.6000.16525" name=GENERATOR>
<STYLE></STYLE>
</HEAD>
<BODY bgColor=#ffffff>
<DIV><FONT face=Arial size=2>Hi Luis,</FONT></DIV>
<DIV><FONT face=Arial size=2> thanks for those clarifications.
I forgot some precisions :)</FONT></DIV>
<DIV>
<DIV><FONT face=Arial size=2>Firstly, I didn't paste the setcenter part of my
code, that was above the pasted code.</FONT></DIV></DIV>
<DIV><FONT face=Arial size=2>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.</FONT></DIV>
<DIV><FONT face=Arial size=2>Each stack contains an overlapping volume,
identified by :</FONT></DIV>
<DIV><FONT face=Arial size=2> -
firstFixedOverlappingBlockImage to lastFixedImage in "top stack"</FONT></DIV>
<DIV><FONT face=Arial size=2> -
firstMovingOverlappingBlockImage to lastMovingImage in "bottom
stack"</FONT></DIV>
<DIV><FONT face=Arial size=2>Thus I can easily geometrically calculate the
center of the overlapping volume in the fixed stack and the same for the moving
stack.</FONT></DIV>
<DIV><FONT face=Arial size=2>Translation is then set to
(centerMovingImage-centerFixedImage) and rotation center set to the
fixedOverlappingVolumeCenter.</FONT></DIV>
<DIV><FONT face=Arial size=2>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..</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial
size=2>********************************************************************************************</FONT></DIV>
<DIV><FONT face=Arial size=2> //
----------------------------------------<BR> // OVERLAPPING BLOCKS CENTER
COMPUTATION<BR> // ----------------------------------------<BR> //
FIXED IMAGE<BR> const SourceImageType::SpacingType fixedSpacing =
pile3Drecto->GetSpacing();</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> double
fixedImageOverlappingBlockCoordinates[3];<BR> fixedImageOverlappingBlockCoordinates[0]
= fixedImageSize[0] / 2.0;<BR> fixedImageOverlappingBlockCoordinates[1] =
fixedImageSize[1] / 2.0;<BR> fixedImageOverlappingBlockCoordinates[2] =
(lastFixedImage+firstFixedOverlappingBlockImage)/ 2.0;</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> TransformType::InputPointType
centerFixedImage;<BR> centerFixedImage[0] = fixedOrigin[0]+ fixedSpacing[0]
* fixedImageOverlappingBlockCoordinates[0];<BR> centerFixedImage[1] =
fixedOrigin[1]+ fixedSpacing[1] *
fixedImageOverlappingBlockCoordinates[1];<BR> centerFixedImage[2] =
fixedOrigin[2]+ fixedSpacing[2] *
fixedImageOverlappingBlockCoordinates[2];</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> // MOVING IMAGE<BR> const
SourceImageType::SpacingType movingSpacing =
pile3Dverso->GetSpacing();</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> double
movingImageOverlappingBlockCoordinates[3];<BR> movingImageOverlappingBlockCoordinates[0]
= (movingImageSize[0] /
2.0)+offsetX;<BR> movingImageOverlappingBlockCoordinates[1] =
(movingImageSize[1] /
2.0)+offsetY;<BR> movingImageOverlappingBlockCoordinates[2] =
((lastMovingImage+firstMovingOverlappingBlockImage)/ 2.0)+offsetZ;</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> TransformType::InputPointType
centerMovingImage;<BR> centerMovingImage[0] = movingOrigin[0]+
movingSpacing[0] *
movingImageOverlappingBlockCoordinates[0];<BR> centerMovingImage[1] =
movingOrigin[1]+ movingSpacing[1] *
movingImageOverlappingBlockCoordinates[1];<BR> centerMovingImage[2] =
movingOrigin[2]+ movingSpacing[2] *
movingImageOverlappingBlockCoordinates[2];<BR> cout << "Centre de
rotation MOVING (coordonnees): [" <<
movingImageOverlappingBlockCoordinates[0] << ";"
<<movingImageOverlappingBlockCoordinates[1] << ";" <<
movingImageOverlappingBlockCoordinates[2] << "]" <<
endl;<BR> cout << "Centre de rotation MOVING (um): [" <<
centerMovingImage[0] << ";" << centerMovingImage[1] << ";"
<< centerMovingImage[2] << "]" << endl;<BR> cout <<
"OFFSET (coordonnees): [" << offsetX << ";" << offsetY
<< ";" << offsetZ << "]" << endl;</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> //
-------------------------------------------------------------------------<BR> //
SET ROTATION CENTER SET AND TRANSLATION VECTOR<BR> //
-------------------------------------------------------------------------<BR> //
set center of fixed image as the rotation center<BR> <FONT
color=#ff0000>transform->SetCenter( centerFixedImage );</FONT><BR> //
translation set as the vector relating the center of the moving image to the
center of the fixed image<BR> <FONT
color=#ff0000>transform->SetTranslation(
centerMovingImage-centerFixedImage);</FONT></FONT></DIV>
<DIV><FONT face=Arial
size=2>********************************************************************************************</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>For the quaternion computation, in fact my real
code is a little bit different than the bad test code I sent.</FONT></DIV>
<DIV><FONT face=Arial size=2>As you said, I constructed the quaternion for each
rotation, multiplied and set the rotation like this:</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial
size=2>*************************************************************</FONT></DIV>
<DIV><FONT face=Arial size=2> typedef TransformType::VersorType
VersorType;<BR> VersorType
rotation;<BR></FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> // rotation d'un angle angleZ selon Oz, dans
les deux cas de flip<BR> double RotZ = angleZ * vnl_math::pi /
180.0;<BR> vnl_quaternion<double> Rz(0, 0,
RotZ);<BR> vnl_quaternion<double> RFlip, RR;</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> if (flipHorizVert) RFlip =
vnl_quaternion<double>(0, vnl_math::pi, 0); // rotation de 180° selon
Oy<BR> else RFlip = vnl_quaternion<double>(vnl_math::pi, 0, 0); //
rotation de 180° selon Ox</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> RR = RFlip *
Rz;<BR> rotation.Set(RR);</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2> transform->SetRotation( rotation
);<BR> registration->SetInitialTransformParameters(
transform->GetParameters() );</FONT></DIV>
<DIV><FONT face=Arial
size=2>*************************************************************</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>All this code gave a quite nice initialization
result:</FONT></DIV>
<DIV><FONT face=Arial size=2>- 3D view: <A
href="http://itk.fete.ch/perso/3DViewInit.jpg">http://itk.fete.ch/perso/3DViewInit.jpg</A></FONT></DIV>
<DIV><FONT face=Arial size=2>- Slice view: <A
href="http://itk.fete.ch/perso/SliceInit001.jpg">http://itk.fete.ch/perso/SliceInit001.jpg</A> and
<A
href="http://itk.fete.ch/perso/SliceInit002.jpg">http://itk.fete.ch/perso/SliceInit002.jpg</A></FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>But as you can see on <A
href="http://itk.fete.ch/perso/SliceInit003.jpg">http://itk.fete.ch/perso/SliceInit003.jpg</A>,
the first (red) image of the moving stack via the initialized transform is half
cut.</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>I'll change my code for your versor code and see if
it happend again.</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>Cheers, nic</FONT></DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2></FONT> </DIV>
<DIV><FONT face=Arial size=2>----- Original Message ----- </FONT>
<DIV><FONT face=Arial size=2>From: "Luis Ibanez" <</FONT><A
href="mailto:luis.ibanez@kitware.com"><FONT face=Arial
size=2>luis.ibanez@kitware.com</FONT></A><FONT face=Arial
size=2>></FONT></DIV>
<DIV><FONT face=Arial size=2>To: "Nic" <</FONT><A
href="mailto:itk@fete.ch"><FONT face=Arial size=2>itk@fete.ch</FONT></A><FONT
face=Arial size=2>></FONT></DIV>
<DIV><FONT face=Arial size=2>Cc: "Insight Users" <</FONT><A
href="mailto:insight-users@itk.org"><FONT face=Arial
size=2>insight-users@itk.org</FONT></A><FONT face=Arial size=2>></FONT></DIV>
<DIV><FONT face=Arial size=2>Sent: Thursday, October 11, 2007 5:44
PM</FONT></DIV>
<DIV><FONT face=Arial size=2>Subject: Re: [Insight-users] Volume not perfectly
rotated</FONT></DIV></DIV>
<DIV><FONT face=Arial><BR><FONT size=2></FONT></FONT></DIV><FONT face=Arial
size=2>> <BR>> Hi Nic,<BR>> <BR>> <BR>> It seems that you have
not initialized the Center of Rotation of the<BR>>
VersorRigid3DTransform.<BR>> <BR>> By default that center is at (0,0,0),
so it is likely that your image<BR>> if being flipped around the corner of
the image, not around the center<BR>> of the image.<BR>> <BR>> You may
want to compute the center of the image in physical coordinates<BR>> and use
that as the center of rotation for the VersorRigid3DTransform.<BR>> <BR>>
You will find examples on how to compute the center of the image in<BR>> the
ITK Software Guide.<BR>> <BR>>
</FONT><A href="http://www.itk.org/ItkSoftwareGuide.pdf"><FONT face=Arial
size=2>http://www.itk.org/ItkSoftwareGuide.pdf</FONT></A><BR><FONT face=Arial
size=2>> <BR>> <BR>> It comes down to<BR>>
<BR>> for( i=0; i<3;
i++)<BR>>
{<BR>> Center[i] = Origin[i] + Size[i] *
Spacing[i] / 2.0<BR>> }<BR>> <BR>>
-----<BR>> <BR>> <BR>> Also, Please aware that *you are not using the
Quaternion correctly*.<BR>> <BR>> <BR>> //: Construct
quaternion from Euler Angles,<BR>> // That is a rotation about
the X axis, followed by Y, followed by<BR>> // the Z axis, using
a fixed reference frame.<BR>> vnl_quaternion(T theta_X, T
theta_Y, T theta_Z);<BR>> <BR>> <BR>> <BR>> That is not the order of
the rotational sequence that you are assuming.<BR>> <BR>> <BR>>
<BR>> It is too bad that so many people have grown used to Euler
angles.<BR>> They are one of the *worst things* that have made their way
in<BR>> computer graphics. They are clumsy and provide a very poor<BR>>
representation of the SO(3) rotational space<BR>>
<BR>> </FONT><A
href="http://en.wikipedia.org/wiki/SO%283%29"><FONT face=Arial
size=2>http://en.wikipedia.org/wiki/SO%283%29</FONT></A><BR><FONT face=Arial
size=2>> <BR>> Euler angles have impoverished the understanding of
generations of<BR>> engineers on the nature of rotational space and crippled
them to<BR>> think in Cartesian schems instead of embracing the natural
properties<BR>> of the rotational space. That results in uncountable pieces
of software<BR>> that require "if" conditions when computing rotational
components.<BR>> <BR>> <BR>> <BR>> ---------<BR>> <BR>> The
correct way of composing rotations in space when using<BR>> Quaternions is to
define a quaternion for each one of the<BR>> intermediate rotations and then
use the Quaternion composition<BR>> method:<BR>>
<BR>> vnl_quaternion::operator*() in
VNL<BR>> <BR>> or<BR>> <BR>>
itkVersor::operator*() in ITK<BR>> <BR>> <BR>> In your
case:<BR>> <BR>> > What I do is simply initialize the
versorRigid3dTransform with a 180°<BR>> > rotation around y-axis and
-94.55 rotation around z-axis.<BR>> <BR>> <BR>> You should do:<BR>>
<BR>> typedef itk::Versor<double>
VersorType;<BR>> typedef VersorType::VectorType AxisType;<BR>>
typedef VersorType::ValueType ValueType;<BR>> <BR>>
AxisType axis1; // Rotation Axis = Y<BR>> axis1[0] = 0.0;<BR>>
axis1[0] = 1.0;<BR>> axis1[0] = 0.0;<BR>> <BR>> ValueType angle1 =
3.141516...<BR>> <BR>> VersorType versor1;<BR>> versor1.Set( axis1,
angle1 ); // 180 degrees around Y<BR>> <BR>> <BR>> AxisType
axis2; // Rotation Axis = Z<BR>> axis2[0] = 0.0;<BR>> axis2[0]
= 0.0;<BR>> axis2[0] = 1.0;<BR>> <BR>> ValueType angle2 = 3.141516 * (
-94.55 ) / 180.0;<BR>> <BR>> VersorType versor2;<BR>> versor2.Set(
axis2, angle2 ); // -94.55 degrees around Z<BR>> <BR>>
<BR>> VersorType versor3 = versor2 *
versor1;<BR>> <BR>> <BR>> The last variable "versor3" will contain the
correct rotation<BR>> expressed in terms of a unit quaternion.<BR>>
<BR>> <BR>> Note that most of the time, when people use
Quaternions,<BR>> what they actually want to use are Versors.<BR>>
<BR>> Quaternions can represent rotation *and scaling* in space,<BR>>
while Versors are limited to rotations. Versors are equivalent<BR>> to
Unit-Quaternions.<BR>> <BR>> <BR>> You may want to read more about
Quaternions and their properties:<BR>> </FONT><A
href="http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation"><FONT
face=Arial
size=2>http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation</FONT></A><BR><FONT
face=Arial size=2>> <BR>> and the ITK Quaternion Tutorials:<BR>>
<BR>> </FONT><A
href="http://www.itk.org/CourseWare/Training/QuaternionsI.pdf"><FONT face=Arial
size=2>http://www.itk.org/CourseWare/Training/QuaternionsI.pdf</FONT></A><BR><FONT
face=Arial size=2>> </FONT><A
href="http://www.itk.org/CourseWare/Training/QuaternionsII.pdf"><FONT face=Arial
size=2>http://www.itk.org/CourseWare/Training/QuaternionsII.pdf</FONT></A><BR><FONT
face=Arial size=2>> <BR>> <BR>> Regards,<BR>>
<BR>> <BR>> Luis<BR>> <BR>>
<BR>> <BR>> ----------------<BR>> Nic wrote:<BR>>> Hi
all,<BR>>> I have a strange behaviour I would like
to correct in my initialisation.<BR>>> I used quaternion for rotation
initialisation like this, as previously<BR>>> adviced :<BR>>>
<BR>>> *******************************************************<BR>>>
double RotX = angleX * vnl_math::pi / 180.0;<BR>>> double RotY = angleY *
vnl_math::pi / 180.0;<BR>>> double RotZ = angleZ * vnl_math::pi /
180.0;<BR>>> vnl_quaternion<double> RR(RotX,RotY,RotZ);<BR>>>
rotation.Set(RR);<BR>>> transform->SetRotation( rotation );<BR>>>
<BR>>> [...]<BR>>> // Checking<BR>>> VersorType
rotationInitiale = transform->GetVersor();<BR>>>
vnl_vector_fixed<double,3> anglesEuleriensInitiaux = <BR>>>
(rotationInitiale.GetVnlQuaternion()).rotation_euler_angles();<BR>>> cout
<< "Angles d'Euler (initiaux): " << <BR>>>
(anglesEuleriensInitiaux*180.0/vnl_math::pi) << endl;<BR>>>
<BR>>>
********************************************************<BR>>>
<BR>>> What I doesn't understand actually is why by flipping around X or Y
axis <BR>>> by 180°, I get an "half image" at the beginning of the
stack<BR>>> Is this behavious linked to the Quaternions ? Is there a
special case <BR>>> for a pi angle ? Is there a way to avoid this
behaviour ?<BR>>> <BR>>> Images:<BR>>> </FONT><A
href="http://itk.fete.ch/perso/Rot_0_0_0.jpg"><FONT face=Arial
size=2>http://itk.fete.ch/perso/Rot_0_0_0.jpg</FONT></A><BR><FONT face=Arial
size=2>>> </FONT><A href="http://itk.fete.ch/perso/Rot_0_180_0.jpg"><FONT
face=Arial size=2>http://itk.fete.ch/perso/Rot_0_180_0.jpg</FONT></A><BR><FONT
face=Arial size=2>>> </FONT><A
href="http://itk.fete.ch/perso/Rot_180_0_0.jpg"><FONT face=Arial
size=2>http://itk.fete.ch/perso/Rot_180_0_0.jpg</FONT></A><BR><FONT face=Arial
size=2>>> <BR>>> <BR>>> <BR>>>
------------------------------------------------------------------------<BR>>>
<BR>>> _______________________________________________<BR>>>
Insight-users mailing list<BR>>> </FONT><A
href="mailto:Insight-users@itk.org"><FONT face=Arial
size=2>Insight-users@itk.org</FONT></A><BR><FONT face=Arial size=2>>>
</FONT><A href="http://www.itk.org/mailman/listinfo/insight-users"><FONT
face=Arial
size=2>http://www.itk.org/mailman/listinfo/insight-users</FONT></A><BR><FONT
face=Arial size=2>></FONT> </BODY></HTML>