[Insight-developers] [Insight-users] Convert Rigid3DTransform To VersorRigid3DTransform
Hans Johnson
hans-johnson at uiowa.edu
Thu Dec 18 12:59:11 EST 2008
Brad,
Thanks for looking into this. I took a less theoretical method in
convincing myself that this was working. I built a test case (
RotationMatrixToVersorTest at line 92 in
Testing/Code/Common/itkVersorTest.cxx ) where I created rotation matrices
to within 1/100000 of a degree around problematic rotations (pi, 0, 90
degrees) then converted those matrices to versors applied them to a point
far away from the origin and measured confirmed that the error in the
separation of rotated points was less than 1e-5/3.
This test was failing with epsilon = 1e-30, but passes with
epsilon=vcl_numeric_limits<double>
The itk::Versor is templated over type, and all the math is done with either
double precision or type T. There are no explict floating point
computations in the conversion from RotationMatrix to Versor.
============
On a separate topic: In the operator== for versor, it seems strange that
for a templated class of the T (floats or doubles) that epsilon is set to
1e-300. That seems too stringent for a floating point comparison.
Additionally, on line 179 the there is a definition of 1.0f, even though the
computation is supposed to be done in double when T=double. There are
other places where double precision is computed with floating point
constants (vcl_sqrt( 2.0f ) instead of vcl_sqrt(2.0) or vnl_math::sqrt2 ).
And finally there are places where division could be replaced by
multiplication (angle/2.0) --> (0.5*angle);
There are three types used in this class, and it seems that they have not
been used consistently: ValueType, T, RealType, double. It seems that all
the references to double should be replaced by RealType in the itkVersor.txx
file, and all references to T should be replaced with ValueType in the
itkVersor.txx file.
Let me know what you think, and I¹ll make this ³Style² changes if you agree
that they will not change the desired outputs.
Hans
On 12/18/08 10:46 AM, "Bradley Lowekamp" <blowekamp at mail.nih.gov> wrote:
> Hello,
>
> It's been a while a while since I used a versor/quaternion much so I am not
> completely sure of some of how this will be used. So this is the bit of code:
>
> const T epsilon = vcl_numeric_limits<T>::epsilon();
> double trace = m(0,0) + m(1,1) + m(2,2) + 1.0;
> if( trace > epsilon)
> {
> const double s = 0.5 / vcl_sqrt(trace);
> m_W = 0.25 / s;
> m_X = (m(2,1) - m(1,2)) * s;
> m_Y = (m(0,2) - m(2,0)) * s;
> m_Z = (m(1,0) - m(0,1)) * s;
> }
> else ..
> Let me see what I can remember from numerical analysis. lf we assume that the
> trace is 2*eps then we will only have about 2 bits of precision, since 1.0 is
> in the computation of trace. This does not look like enough precision. The
> other methods may be better before this point. Perhaps eps should be
> multiplied by the number of bits needed of the computation. Hrmm.. I am
> uncertain, but I think some of that computation is only occurring as floats
> too.
>
> Just some thought for consideration :)
>
> Brad
>
> On Dec 17, 2008, at 11:41 PM, Hans Johnson wrote:
>
>> All,
>>
>> I wanted to follow up a bit on what I've determined with regards to
>> converting between itkVersor and Rotation Matrices (and back again). The
>> itk implmentation had the threshold for using alternate computations to
>> compute the conversions near the singularities that arise near roations of
>> pi as 1e-30, but at that very small value, the numerical stability of the
>> standard calculation was not sufficient to compute proper matricies. Upon
>> reading
>> http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternio
>> nToEuler/ I determined that the threshold needed to be changed.
>>
>> A bug has been submitted to http://public.kitware.com/Bug/view.php?id=8312
>>
>> I've changed this threshold to vcl_numeric_limits<double>::epsilon(), and
>> wrote a test to confirm that the proper desired behavior was occuring.
>>
>> /cvsroot/Insight/Insight/Code/Common/itkVersor.txx,v <--
>> Code/Common/itkVersor.txx
>> new revision: 1.29; previous revision: 1.28
>> /cvsroot/Insight/Insight/Testing/Code/Common/itkVersorTest.cxx,v <--
>> Testing/Code/Common/itkVersorTest.cxx
>> new revision: 1.20; previous revision: 1.19
>>
>> ============================================
>> I believe that this change (and the test case that demonstrates it is stable
>> and reliable) allows itkVersorRigid3DTransform.h SetRotationMatrix function
>> to be declared as a safe way to initialize a the Versor.
>>
>> Any comments would be greatly appreciated.
>>
>> Thanks,
>> Hans
>>
>>
>>
>> On 12/17/08 7:32 AM, "Hans Johnson" <hans-johnson at uiowa.edu> wrote:
>>
>>> Bill,
>>>
>>> Thanks. I have that page printed on my desk and have been referencing it.
>>>
>>> The code I have supplied in my previous message is producing a
>>> VersorRigid3DTransformation that produces different resampling results than
>>> the original Rigid3DTransform.
>>>
>>> ========= From itkVersor.h ============
>>> 205 /** Set the versor using an orthogonal matrix.
>>> 206 * Based on code from:
>>> 207 * http://www.euclideanspace.com/maths/geometry/rotations/
>>> 208 * conversions/matrixToQuaternion/index.htm
>>> 209 */
>>> 210 void Set( const MatrixType & m );
>>>
>>> ======== From itkVersorRigid3DTransform.h ================
>>> 116
>>> 117 /** This method must be made protected here because it is not a safe
>>> way of
>>> 118 * initializing the Versor */
>>> 119 virtual void SetRotationMatrix(const MatrixType & matrix)
>>> 120 { this->Superclass::SetRotationMatrix( matrix ); }
>>>
>>>
>>> The comment in line 117 of itkVersorRigid3DTransform.h is concerning to me.
>>> I would think that this would be a safe way to initialize the rotation part
>>> of the itkVersor3DTransform, as long as the input matrix was tested to be
>>> orthogonal.
>>>
>>> ======== From itkVersorRigid3DTransform.h ================
>>> 70 // Check if input matrix is orthogonal to within tolerance
>>> 71 template<class TScalarType>
>>> 72 bool
>>> 73 Rigid3DTransform<TScalarType>
>>> 74 ::MatrixIsOrthogonal(
>>> 75 const MatrixType & matrix,
>>> 76 double tolerance )
>>> 77 {
>>> 78 typename MatrixType::InternalMatrixType test =
>>> 79 matrix.GetVnlMatrix() * matrix.GetTranspose();
>>> 80
>>> 81 if( !test.is_identity( tolerance ) )
>>> 82 {
>>> 83 return false;
>>> 84 }
>>> 85
>>> 86 return true;
>>> 87 }
>>>
>>>
>>> The itkVersorRigid3DTransform has 6 total parameters. The last three are
>>> the translation, and the first 3 are the rotation. The versor is a unit
>>> quaternion, so one of the 4 quaternion parameters can be computed from the
>>> other 3. I am not sure if I am placing the proper 3 parameters in the first
>>> three locations of the itk::VersorRigid3DTransform.
>>>
>>> I would like to see the following code work, or an explanation of why this
>>> will not work:
>>>
>>> itk::VersorRigid3DTransform<double>::Pointer
>>> ConvertToVersorRigid3D( itk::Rigid3DTransform<double>::Pointer RT)
>>> {
>>> if(! RT->MatrixIsOrthogonal(RT->GetRotationMatrix()))
>>> {
>>> throw exception...
>>> }
>>> itk::VersorRigid3DTransform<double>::Pointer returnVersorRigid3D=
>>> itk::VersorRigid3DTransform<double>::New();
>>> returnVersorRigid3D->SetFixedParameters(RT->GetFixedParameters());
>>> returnVersorRigid3D->SetRotationMatrix( RT->GetRotationMatrix() );
>>> returnVersorRigid3D->SetTranslation( RT->GetTranslation() );
>>> return returnVersorRigid3D;
>>> }
>>>
>>>
>>>
>>> Thanks,
>>> Hans
>>>
>>>
>>>
>>>
>>> On 12/16/08 11:02 PM, "Bill Lorensen" <bill.lorensen at gmail.com> wrote:
>>>
>>>> Hans,
>>>>
>>>> Here's a wikipedia entry:
>>>>
http://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angle>>>>
s
>>>>
>>>> Bill
>>>>
>>>> On Tue, Dec 16, 2008 at 10:20 PM, Hans Johnson <hans-johnson at uiowa.edu>
>>>> wrote:
>>>>> Hello All,
>>>>>
>>>>> I have a program that generates rho, phi, theta rotation euler angles, and
>>>>> I
>>>>> can easily generate the rotation matrix for creating the Rigid3DTransform.
>>>>>
>>>>> I now need to convert Rigid3DTransform to VersorRigid3DTransform so that I
>>>>> can initialize a mutual information registration with it.
>>>>>
>>>>> Could anyone provide some advice on how to perform this conversion?
>>>>>
>>>>> Thanks,
>>>>> Hans
>>>>>
>>>>> PS: The following code does not work.
>>>>> ===========================
>>>>> typedef itk::Rigid3DTransform<double> RigidTransformType;
>>>>> typedef itk::VersorRigid3DTransform<double> VersorTransformType;
>>>>>
>>>>> VersorTransformType::Pointer
>>>>> ConvertToVersorRigid3D(RigidTransformType::Pointer RT)
>>>>> {
>>>>> VersorTransformType::Pointer VT=VersorTransformType::New();
>>>>> VT->SetFixedParameters(RT->GetFixedParameters());
>>>>>
>>>>> itk::Matrix<double,3,3> R=RT->GetRotationMatrix();
>>>>> RigidTransformType::TranslationType T=RT->GetTranslation();
>>>>>
>>>>> VersorTransformType::ParametersType p;
>>>>> p.SetSize(6);
>>>>> itk::Versor<double> v;
>>>>> v.Set(R);
>>>>> //Get the first 3 elements of the versor;
>>>>> p[0]=v.GetRight()[0];
>>>>> p[1]=v.GetRight()[1];
>>>>> p[2]=v.GetRight()[2];
>>>>> p[3]=T[0];
>>>>> p[4]=T[1];
>>>>> p[5]=T[2];
>>>>> VT->SetParameters(p);
>>>>> return VT;
>>>>> }
>>>>>
>>>>>
>>>>> Notice: This UI Health Care e-mail (including attachments) is covered by
>>>>> the
>>>>> Electronic Communications Privacy Act, 18 U.S.C. 2510-2521, is
>>>>> confidential
>>>>> and may be legally privileged. If you are not the intended recipient, you
>>>>> are hereby notified that any retention, dissemination, distribution, or
>>>>> copying of this communication is strictly prohibited. Please reply to the
>>>>> sender that you have received the message in error, then delete it. Thank
>>>>> you.
>>>>>
>>>>> _______________________________________________
>>>>> Insight-users mailing list
>>>>> Insight-users at itk.org
>>>>> http://www.itk.org/mailman/listinfo/insight-users
>>>>>
>>>>>
>>>
>>> _______________________________________________
>>> Insight-users mailing list
>>> Insight-users at itk.org
>>> http://www.itk.org/mailman/listinfo/insight-users
>>
>> _______________________________________________
>> Insight-users mailing list
>> Insight-users at itk.org
>> http://www.itk.org/mailman/listinfo/insight-users
>
>
> ========================================================
>
> Bradley Lowekamp
>
> Lockheed Martin Contractor for
>
> Office of High Performance Computing and Communications
>
> National Library of Medicine
>
> blowekamp at mail.nih.gov
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20081218/1310b033/attachment.htm>
More information about the Insight-developers
mailing list