[Insight-users] to understand quaternion versor and matrix in VersorRigid3DTransform
Serena Fabbri
fabbri at u.washington.edu
Fri Mar 20 18:10:21 EDT 2009
Hi Luis, hi All,
thank for your exhaustive explanation.
I think wiki-link got mixed up me.
I thought that versor was rotating around Z as I set Z like rotating axes in my code and if i see the m33 element of rotation matrix is very close to 1...so i guessed it.
but now i know that the rotation matrix is the one in the PDF document of the Quaternions tutorial.
Thanks a lot!!!
Serena
------------------------------
>
> Message: 6
> Date: Fri, 20 Mar 2009 15:03:29 -0400
> From: Luis Ibanez <luis.ibanez at kitware.com>
> Subject: Re: [Insight-users] to understand quaternion versor and
> matrix in VersorRigid3DTransform
> To: Serena Fabbri <fabbri at u.washington.edu>
> Cc: itk <insight-users at itk.org>
> Message-ID:
> <f7abd23c0903201203s536f4f30m359b39e2e4e9bb83 at mail.gmail.com>
> Content-Type: text/plain; charset=ISO-8859-1
>
> Hi Serena,
>
>
> Thanks for your detailed question.
>
>
> 1) The Versor is a unit Quaternion.
>
> The reason why we can represent it with only three numbers
> (instead of the four in the quaternion) is that the fourth number
> is already determined by the requirement of the Versor norm
> to be == 1.0.
>
> This means that if you have a Versor of components x,y,z,
> the fourth component w is fully defined by:
>
> w^2 = 1.0 - x^2 - y^2 - z^2
>
> It is therefore, perfectly normal to find that
>
> x^2 + y^2 + z^2 < 1.0
>
> what should never happen is:
>
> x^2 + y^2 + z^2 > 1.0
>
>
> 3) How to interpret a Versor ?
>
> Here is a practical way to do it:
>
> a) compute K = sqrt( x^2 + y^2 + z^2 )
>
> K is equal to: sin( theta / 2.0 )
>
> where theta is the angle of rotation.
> in other words:
>
> theta = 2.0 * asin( K )
>
> b) take the x,y,z components and divide
> them by K. That will give you the unit
> vector parallel to the rotation axis.
>
> Axis = ( x/K, y/K, z/K )
>
> In your particular case,
>
>> versor X = -0.00431463
>> versor Y = -0.00118787
>> versor Z = 0.00404881
>>
>
> K = 6.0349e-03
>
> therefore
>
> theta = 2.0 * math.asin( 6.0349e-03 ) [radians]
>
> theta = 0.01206 radians = 0.69 degrees
>
>
> Your Versor is not rotating around the Z axis.
> such rotation will have the x and y components
> equal to zero.
>
> Your axis of rotation has components:
>
> Ax = -0.00431463 / K
> Ay = -0.00118787 / K
> Az = 0.00404881 / K
>
>
> 3) The correct rotation matrix is the one in the PDF
> document of the Quaternions tutorial.
>
> You are correct in that *IF* your versor were
> rotating around the Z axis, then the rotation matrix
> will have the form:
>
> C -S 0
> S C 0
> 0 0 1
>
> but, helas,... your versor is not rotating around Z.
>
>
> So, I guess the background question is: Why lead you
> to think that the versor was rotating around Z ?
>
> or is this a requirement of the registration problem
> that you are trying to solve ?
>
>
> Please let us know,
>
> Thanks
>
>
> Luis
>
>
>
> -----------------------------------------------------
> On Fri, Mar 20, 2009 at 2:45 PM, Serena Fabbri <fabbri at u.washington.edu> wrote:
>> Hi All,
>>
>> I have some questions about VersorRigid3DTransform, in particular about the
>> versor and the matrix in the output.
>> I am registering T2 and T1 image.
>>
>> FIX IMAGE: T2 size=(392, 512, 160)
>> spacing=(0.5, 0.5, 1)mm
>>
>> MOVING IMAGE: T1
>> size=(176, 256, 160)
>> spacing=(1, 1, 1)mm
>>
>> I use
>> VersorRigid3DTransform
>> itkMattesMutualInformationImageToImageMetric
>> itkVersorRigid3DTransformOptimizer
>>
>>
>>
>> I find this output:
>>
>>
>> Result =
>> ?versor X ? ? ?= -0.00431463
>> ?versor Y ? ? ?= -0.00118787
>> ?versor Z ? ? ?= 0.00404881
>>
>> ?Translation X = -9.59506
>> ?Translation Y = -14.2873
>> ?Translation Z = 1.52361
>>
>> ?Iterations ? ?= 81
>> ?Metric value ?= -0.204678
>>
>>
>> Matrix = 0.999964 -0.00808722 -0.00241064
>> 0.00810772 0.99993 0.00861949
>> 0.00234077 -0.00863873 0.99996
>>
>> Offset = [-8.36679, -15.7562, 2.40158]
>>
>>
>> The versor is a three first component of Quaternion.
>>
>> In http://en.wikipedia.org/wiki/Versor
>> quaternion is:
>> q = cos(T/2) + u sin(T/2)
>>
>> where u is a unit vector.
>>
>> I have calculated the norm of my versor and it is ||versor||= 6.0349e-03
>> so it is not 1 !!!
>>
>> 1) so...what is the correct interpretation of versor in ITK?
>>
>> 2) In the Luis Ibanez email (Apr 7 2007) , he writes that the quaternion is:
>>
>> ? ? ? q0 = Ax * sin( T / 2 )
>> ? ? ? q1 = Ay * sin( T / 2 )
>> ? ? ? q2 = Az * sin( T / 2 )
>> ? ? ? q3 = ? ? ?cos( T / 2 )
>>
>> where (Ax,Ay,Az) are the components of the axis
>> of rotation and T is the angle of rotation.
>> Versors are composed of only the components
>>
>> ? ? ? ? ?q0,q1,q2
>>
>>
>> I am rotating around z axes so my A components are 0 0 1 and I get
>> q2 = sinT/2 and I guess q2 = versorZ so my final rotation angle is 0.464
>> ?degrees Is this the rotation angle?
>>
>> 3)Is the rotation matrix ?the matrix in slide n. 61 of
>> ?http://www.itk.org/CourseWare/Training/QuaternionsI.pdf ? ?
>>
>> or it is the follow matrix:
>> (to be cleare I am rotating around z axes)
>>
>> M= ?cosT ?sinT 0
>> ? -sinT ?cosT 0
>> ? ? ?0 ? ? 0 ?1
>>
>> because If you read the ITKSoftwareGuide pag.372, to calculate the final
>> rotation angle, it is sufficiently to do arcsin(sinT).
>>
>>
>> Thanks a lot for all explanation.
>>
>>
>> Serena Fabbri
>>
>>
>>
>>
>>
>> _____________________________________
>> Powered by www.kitware.com
>>
>> Visit other Kitware open-source projects at
>> http://www.kitware.com/opensource/opensource.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
>>
>
>
> ------------------------------
>
> Message: 7
> Date: Fri, 20 Mar 2009 15:31:50 -0400
> From: Luis Ibanez <luis.ibanez at kitware.com>
> Subject: [Insight-users] YouTube - Authors at Google: Donald Knuth
> To: Insight Users <insight-users at itk.org>
> Message-ID: <49C3EF26.7020506 at kitware.com>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
>
> Donald Knuth talk at Googleplex:
>
> http://www.youtube.com/watch?v=JPpk-1btGZk&feature=sdig&et=1237572947.75
>
> http://en.wikipedia.org/wiki/Donald_Knuth
>
> Author of
>
> "The Art of Computing Programming"
>
> and...
>
> TeX
>
>
>
>
>
> ------------------------------
>
> Message: 8
> Date: Fri, 20 Mar 2009 15:54:52 -0400
> From: Michael Jackson <mike.jackson at bluequartz.net>
> Subject: [Insight-users] Questions about Registering 2 image with
> small overlap
> To: insight-users at itk.org
> Message-ID: <B22BBE97-AF86-4F23-A11F-849FA37B6B8C at bluequartz.net>
> Content-Type: text/plain; charset=US-ASCII; format=flowed; delsp=yes
>
> I'm _still_ new to the whole ITK/Image Registration thing. I have been
> over the ITK manual a few times (selected sections that I thought was
> relevant) but I am having a tough time getting a registration to work
> out. Here is what I have.
>
> 2 Images. Both are 1292 pixels wide by 968 pixels high. The scaling of
> the image is 0.207987 microns/pixel. The origins (upper left, not the
> it probably matters) are:
>
> Image 0: 47113.2, 48448.5 (microns)
> Image 1: 47369.2, 48448.8 (microns).
>
> So basically the edges overlap by a little bit (about 20 pixels or
> so). The edges are not quite lined up properly with those settings
> (which came from the instrument the image was captured on) so I am
> trying to run the images through ITK in order to align them better.
>
> I am getting the exception: Description: itk::ERROR:
> MeanSquaresImageToImageMetric(0x569100): Too many samples map outside
> moving image buffer: 13840 / 1722508
>
> which I _think_ is telling me that something is "off" in where the
> images are being placed in physical space.
>
> I am not even sure if I am using the proper Image Metric for this type
> of problem. The images are grayscale in nature.
>
> Here is some code for your perusal. Any pointers on what I might be
> doing wrong would surely be appreciated.
>
> Thanks
> Mike Jackson
>
>
> const unsigned int Dimension = 2;
> typedef unsigned char PixelType;
> typedef itk::Image<PixelType, Dimension> FixedImageType;
> typedef itk::Image<PixelType, Dimension> MovingImageType;
> typedef itk::TranslationTransform<double, Dimension> TransformType;
> typedef itk::RegularStepGradientDescentOptimizer
> OptimizerType;
> typedef itk::MeanSquaresImageToImageMetric
> <FixedImageType, MovingImageType> MetricType;
> typedef itk::LinearInterpolateImageFunction
> <MovingImageType, double> InterpolatorType;
> typedef itk::ImageRegistrationMethod
> <FixedImageType, MovingImageType> RegistrationType;
> typedef itk::ImageFileReader< FixedImageType > ReaderType;
> typedef itk::ImageFileWriter< FixedImageType > WriterType;
> typedef itk::Point< double, FixedImageType::ImageDimension > PointType;
>
>
> FixedImageType::PointType fixedOrigin;
> fixedOrigin[0] = 47113.2;
> fixedOrigin[1] = 48448.5;
>
> MovingImageType::PointType movingOrigin;
> movingOrigin[0] = 47369.2;
> movingOrigin[1] = 48447.8;
>
> fixedImageReader->GetOutput()->SetOrigin(fixedOrigin);
> movingImageReader->GetOutput()->SetOrigin(movingOrigin);
>
> //--- Set the Spacing of the images
> FixedImageType::SpacingType fixedSpacing;
> fixedSpacing[0] = 0.207987;
> fixedSpacing[1] = 0.207987;
> MovingImageType::SpacingType movingSpacing;
> movingSpacing[0] = 0.207987;
> movingSpacing[1] = 0.207987;
>
> fixedImageReader->GetOutput()->SetSpacing(fixedSpacing);
> movingImageReader->GetOutput()->SetSpacing(movingSpacing);
>
> registration->SetFixedImage(fixedImageReader->GetOutput());
> registration->SetMovingImage(movingImageReader->GetOutput());
>
> registration->SetFixedImageRegion(
> fixedImageReader->GetOutput()->GetBufferedRegion());
>
> typedef RegistrationType::ParametersType ParametersType;
> ParametersType initialParameters(transform->GetNumberOfParameters());
> initialParameters[0] = INITIAL_TRANS_X; // Initial offset in mm along X
> initialParameters[1] = INITIAL_TRANS_Y; // Initial offset in mm along
> Y registration->SetInitialTransformParameters(initialParameters);
>
>
>
>
> ------------------------------
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
>
> End of Insight-users Digest, Vol 59, Issue 70
> *********************************************
>
More information about the Insight-users
mailing list